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CHAPTER  I 
INTRODUCTION 


1.1  Motivation,  Scope  and  Organization 
1.1.1  Main  Objective 

The  main  objective  of  this  research  is  to  develop  and  demonstrate 
a  systematic  and  unified  theory  and  computational  method  for  the 
design  of  large  scale  constrained  dynamic  mechanisms  and  machines.  The 
key  to  meeting  this  objective  lies  in  judicious  selection  of  the  most 
suitable  methods  from  the  following  branches  of  mathematics  and  mechan¬ 
ics  : 


(a)  Methods  of  Optimization 

(b)  Rigid  Body  Mechanics 

(c)  Numerical  Integration  Methods 

(d)  Matrix  Manipulation  Methods. 

The  objective  can  be  illustrated  by  the  following  Venn  diagram: 


Figure  1.1  Venn  Diagram  for  the  Main  Objective 
of  the  Research. 
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The  arrows  on  the  diagram  may  be  viewed  as  the  process  of  selecting 
from  among  a  large  number  of  alternatives  in  each  area  a  method  that  is 
compatible  with  those  of  the  other  areas,  together  providing  a  qualita- 
tively  new  design  capability.  Such  a  comprehensive  treatment  has  never 
been  attempted,  so  utmost  possible  care  has  been  taken  in  the  theoreti¬ 
cal  investigation  to  implement  the  most  suitable  computer  algorithms 
associated  with  various  branches  of  mechanics  and  mathematics  noted 
above.  Brief  discussions  of  the  above  branches,  together  with  the 
indication  of  the  methods  selected,  are  given  below. 

1.1.2  The  Notion  of  Optimal  Engineering  Design 

The  job  of  "optimal  engineering  design"  is  to  develop  the  best 
possible  system  for  the  given  application,  consistent  with  resources 
allocated  to  the  development  phase.  Although  the  notion  of  optimiza¬ 
tion  is  inherent  in  the  design  process,  optimization  as  a  formalized 
approach  to  engineering  design  is  a  relatively  new  concept.  It  is  in 
the  judicious  selection  of  a  "quantitative  measure  of  performance"  of 
the  system  and  in  quantifying  performance  constraints  that  an  optimal 
design  is  distinguished  from  a  conventional  design. 

After  quantification  of  the  notions  of  a  design  process  comes  the 
role  of  selection  of  numerical  methods  for  its  solution,  using  the 
modern  high-speed  digital  computer.  Great  strides  have  been  made  with 
digital  computers  in  the  past  two  decades,  to  allow  for  numerical 
analysis  as  a  test  of  an  idea  or  concept,  rather  than  previous  cut-and- 
try  techniques.  In  this  report,  advantage  is  taken  of  advances  in 
computer-aided  design  techniques  [1]. 
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The  field  of  optimal  design  of  constrained  dynamic  systems  is 
of  growing  interest  and  importance.  In  the  realm  of  dynamic  mechanisms, 
both  continuous  (smooth)  and  intermittent  (discontinuous)  motion 
occurs.  For  many  mechanisms,  the  logical  sequence  of  events  in  inter¬ 
mittent  motion  is  known  in  advance  and  the  total  period  of  motion  under 
consideration  can  be  divided  into  intervals  of  continuous  motion.  Also, 
some  intermittent  motion  can  be  reduced  to  continuous  motion  by  the 
introduction  of  artificial  spring-damper  systems  (see  Chapter  VI). 

Thus,  derivation  of  a  unified  technique  for  optimal  design  of  dynamic 
systems  with  continuous  motion  is  of  major  importance. 

1.1.3  Methods  of  Optimization 

The  fundamental  problem  of  infinite  dimensional  optimal  design 
can  be  described  by  the  problem  of  Bolza  and  its  extended  version  that 
accounts  for  inequality  functional  constraints  (see  [1]).  A  corre¬ 
sponding  problem  for  constrained  dynamic  systems  with  continuous  motion 
is  formulated  in  Chapter  III.  There  are  many  indirect  methods  [1] 
based  on  a  powerful  Functional  Analysis  theorem  of  Liusternik  and 
Sobolev  [69]  for  the  solution  of  such  problems.  But  generally  they 
pose  serious  computational  hazards.  In  this  dissertation  a  direct 
numerical  method  of  solution  is  adopted. 

The  basic  idea  of  the  direct  method  of  solving  optimal  design 
problems  is  to  first  construct  an  initial  estimate  of  the  solution  and 
then  to  find  small  changes  in  the  design  parameters  such  that  the  modi¬ 
fied  design  forms  an  improved  estimate  of  the  optimum,  in  some  sense. 
Before  design  improvements  can  be  determined,  analysis  of  their  effect 
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on  the  problem  must  be  performed.  This  forms  the  sensitivity  analysis 
part  of  the  design  process  and  is  of  extreme  importance. 

In  the  approach  used  here,  design  sensitivity  analysis  uses 
state-space  methods  (see  [24])  in  which  the  state  variables  and  design 
parameters  are  first  treated  as  independent  variables  and  the  elimina¬ 
tion  of  the  variations  of  the  state  variables  is  performed  through  use 
of  adjoint  differential  and  algebraic  equations  (see  Chapter  III  and 
references  [1,70]). 

1.1.4  Rigid  Body  Mechanics 

There  are  two  general  approaches  to  the  subject  of  rigid  body 
mechanics.  They  can  be  called  the  "Vectorial  or  Newtonian  approach” 
and  the  "Analytical  approach".  Vectorial  dynamics  is  based  on  a  direct 
application  of  Newton's  laws  of  motion  and  concentrates  on  the  forces 
and  motion  associated  with  individual  parts  of  the  system,  whereas 
analytical  dynamics  is  concerned  with  the  system  as  a  whole  and  uses 
descriptive  scalar  functions  such  as  kinetic  and  potential  energies. 

The  most  direct  analytical  approach  is  the  well-known  Lagrangian  for¬ 
mulation.  For  details,  references  [66,67,68]  are  recommended. 

In  most  treatments  of  optimal  design  of  dynamic  systems  in  the 
literature,  the  number  of  generalized  coordinates  of  a  system  is  taken 
equal  to  the  number  of  degrees  of  freedom  of  the  whole  system,  so  that 
with  the  application  of  Newton's  laws,  the  number  of  first  order  ordi¬ 
nary  differential  equations  for  the  state  is  equal  to  twice  the  number 
of  physical  degrees  of  freedom.  In  general,  a  dynamic  mechanism  may 
be  very  complicated  and  direct  application  of  Newton's  laws  may  result 
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in  a  highly  nonlinear  set  of  differential  equations.  Moreover,  the 
eigenvalues  for  such  a  nonlinear  system  may  vary  over  extreme  ranges 
during  the  process  of  numerical  integration.  Therefore,  there  remains 
a  danger  of  the  numerical  integration  problem  turning  "stiff"  [2,3,4, 
5,6]. 

Orlandea  [7]  and  Orlandea,  Chace,  and  Calahan  [8]  have  developed 
a  node-analogous  sparsity-oriented  approach  to  the  dynamic  analysis  of 
mechanical  systems  using  the  Lagrangian  formulation  of  rigid  body 
mechanics.  By  applying  SPARSE  MATRIX  [9,64,65]  and  STIFF  [10]  inte¬ 
gration  algorithms,  large  sets  of  sparse  linear  equations  can  be  effi¬ 
ciently  solved  with  moderate  expenditure  of  CPU  time  and  the  numerical 
instability  [2,3,6,11]  associated  with  widely  split  eigenvalues  at  any 
stage  can  be  avoided.  Both  sparse  matrix  techniques  and  stiff  inte¬ 
gration  algorithms  are  in  the  constant  process  of  modification  and 
refinement.  These  topics  are  treated  briefly  in  Chapter  II  of  this 
report  (also  see  Sections  1.1.5  and  1.1.6). 

These  algorithms  have  been  implemented  by  Orlandea  [7]  to  generate 
a  computer  program  "ADAMS"  (Automatic  Dynamic  Analysis  of  Mechanical 
Systems).  This  program  was  developed  for  efficient  simulation  of 
dynamic  mechanical  systems  (e.g.,  vehicles  and  machinery)  using  methods 
of  numerical  analysis  developed  for  electrical  circuits  [10,12,13]. 
Orlandea' s  work  has  demonstrated  that  the  sparse  matrix  formulation  and 
numerical  methods  involving  stiff  numerical  integration  techniques  can 
be  effectively  used  for  simulation  of  large  three-dimensional  mechanical 
systems.  The  advantages  and  disadvantages  of  these  can  be  stated  as 


.  follows: 
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Advantages: 

1)  No  topological  preprocessing  is  necessary  to  establish  a  set  of 
independent  variables;  equations  can  be  developed  directly  from 
the  connection  data,  component  by  component. 

2)  High  sparsity  of  the  "Jacobian  Matrix"  is  used  in  a  stiff  integra¬ 
tion  algorithm  (see  Chapter  II) . 

3)  All  angular  and  displacement  variables  are  retained  as  solution 
variables;  none  are  eliminated  in  the  interest  of  producing  a 
reduced  set  of  equations  with  fewer  variables.  In  this  way  the 
total  number  of  matrix  operations  and  the  number  of  operations 
for  eliminating  variables  from  nonlinear  equations  are  kept  at  a 
minimum. 

A)  All  joint  reaction  forces  are  determined  directly  in  the  solution 
and  therefore  the  formulation  is  compatible  with  current  methods 
of  continuum  mechanics  for  internal  stress-analysis. 

5)  Frictional  effects  in  joints  are  routinely  handled. 

Disadvantages : 

1)  Nodal  formulation  results  in  more  equations  than  loop  formulation. 

2)  Some  time  may  be  wasted  in  solving  for  variables  of  little  interest 
to  the  designer. 

ADAMS,  IMP  (Integrated  Mechanism  Program)  [15,16],  and  DRAM 
(Dynamic  Response  of  Articulated  Machineries)  [17,18]  are  the  three 
main  computer  programs  at  present  for  dynamic  simulation  of  mechanical 
systems.  Both  DRAM  and  IMP  use  relative  coordinate  systems  for  parts 
(or  bodies  or  nodes)  of  a  mechanism  and  consider  independent  loop 
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equations  for  the  determination  of  constraint  equations.  These  are 
complicated  mathematical  relations  involving  numerous  matrix  inversions. 
ADAMS,  on  the  other  hand,  deals  with  a  global  inertial  reference  frame 
for  all  bodies  in  the  system.  It  will  be  evident  from  the  formulation 
and  analytical  procedure  in  Chapters  II  and  III  that  the  ADAMS  program 
is  extremely  suitable  for  optimal  design  investigations.  All  the  above 
mentioned  programs,  however,  use  Lagrange's  equations  of  motion  in  the 
analysis. 

There  are  other  dynamic  analysis  programs,  namely,  IMP-UM  (IMP  at 
the  University  of  Michigan) ,  MEDUSA  (Machine  Dynamic  Universal  System 
Analyzer)  [19],  VECNET  (Vector  Network)  [20],  and  DYMAC  (Dynamics  of 
Machinery)  [21].  They  deal  with  different  methods  of  rigid  body 
dynamics.  A  more  detailed  overview  of  all  the  programs  can  be  obtained 
in  reference  [22]. 

It  should  be  noted  that,  theoretically,  any  of  the  dynamic  anal¬ 
ysis  programs  discussed  in  the  foregoing  can  be  used  to  predict  system 
dynamics.  An  apparent  advantage  inherent  in  all  the  programs  except 

ADAMS  lies  in  the  fact  that  they  deal  with  a  minimum  number  of  inde¬ 

pendent  (generalized)  coordinates  and  thus  involve  fewer  equations 
than  does  the  ADAMS  method.  The  disadvantages  of  these  methods  are 

much  more  serious,  however.  The  reduced  set  of  equations  with  fewer 

variables  is  highly  nonlinear.  Consequently,  the  eigenvalues  for  such 
nonlinear  systems  may  vary  quite  unpredictably  during  the  integration 
process.  Another  major  disadvantage  lies  in  the  fact  that  they  do  not 
evaluate  the  reaction  forces  simultaneously  with  the  solution  variables. 


In  a  reasonable  optimal  design  formulation,  one  generally  places  bounds 
on  reaction  forces  at  joints  (see  Chapters  III  and  VI).  All  these 
problem  areas  are  readily  handled  by  the  ADAMS  modeling  method,  which 
makes  the  ADAMS  method  extremely  suitable  for  optimal  design  investi¬ 
gations. 

Wehage  [14]  has  written  a  program  for  planar  systems  using  con¬ 
strained  system  formulation  and  sparse  matrix  methods.  This  program, 
after  several  modifications,  has  been  extended  and  implemented  as  the 
analysis  module  of  the  computer  program  "DADS"  (Dynamic  Analysis  and 
Design  System)  in  this  dissertation  (see  Chapters  II,  V). 

1.1.5  Numerical  Integration  Methods 

The  solution  of  linear  dynamic  state  equations  can  always  be 
expressed  in  analytic  forms.  But  when  the  equations  are  nonlinear, 
the  analytic  forms  of  solutions  seldom  exist  and  one  is  compelled  to 
resort  to  graphical  [42]  or  numerical  methods.  The  most  serious  short¬ 
coming  of  graphical  methods  lies  in  their  inapplicability  for  higher 
order,  nonlinear  mechanical  systems.  On  the  other  hand,  numerical 
methods  are  valid  and  appropriate  for  nonlinear  systems  of  any  order. 

In  view  of  their  greater  generality  and  ease  of  implementation  on  a 
digital  computer,  numerical  methods  are  widely  applied.  In  the  follow¬ 
ing,  most  discussions  will  be  confined  to  nonlinear  systems.  For 
detailed  analytic  theory  of  numerical  methods,  references  [2,3,4,6,11, 
43,44,45]  are  recommended. 

In  the  Lagrangian  formulation  of  the  equations  of  motion  of  a 
constrained  mechanical  system,  when  all  the  generalized  coordinates 
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and  Lagrange  multipliers  are  taken  to  be  independent  [66,67],  nonlinear 
algebraic  equations  enter  into  the  system  of  equations  (see  Section  2.1). 
As  put  forward  by  Gear  [46],  such  a  system  can  be  written  as: 

Fci.z.t)  =  o  d.i) 

where  y  is  the  solution  variable  vector,  t  is  the  time  parameter,  and 
y  is  the  vector  of  time  derivatives.  A  member  of  Eq.  (1.1)  may  be 
either  a  differential  or  an  algebraic  equation,  according  a  8F/ 8y  is 
nonzero  or  zero.  Such  a  system  of  equations  is  called  a  simultaneous 
system  of  Differential  and  Algebraic  Equations (DAE' s) .  To  the  system 
of  equations  (1.1)  one  must  add  the  appropriate  initial  conditions  of 
the  form 


i(o)  =  iG  (i-2) 

where  y  is  a  vector  of  the  subset  of  solution  variables  having  time 
derivatives  in  the  system  of  Eq.  (1.1). 

Although  there  exist  many  algorithms  for  the  solution  of  initial 
value  problems,  most  of  them  are  based  on  two  basic  approaches: 

(1)  the  Taylor  series  expansion  approach  and  (2)  polynomial  approxi¬ 
mation  approach.  Algorithms  based  on  the  first  approach  are  generally 
called  Runge-Kutta  algorithms  and  those  based  on  the  second  one  are 
usually  called  numerical  integration  algorithms  (see  references  [2,6]). 
The  former  are  single-step  algorithms,  whereas  the  latter  are  multistep 
algorithms  that  use  information  from  previous  time  steps.  Adams- 
Bashforth  and  Adams-Moulton  algorithms  are  examples  of  the  second 
category  [2,6]. 
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In  most  of  the  treatments  found  in  the  literature,  these  algo¬ 
rithms  deal  with  the  system  of  equations  of  the  form 


i  -  f(x,t) 

z<0)  -  z0 


(1.3) 


which  do  not  include  algebraic  equations.  Gear  [2,46,48]  appears  to 
be  the  first  person  to  present  a  multistep  algorithm  that  deals  with 
a  simultaneous  system  of  differential  and  algebraic  equations,  together 
with  the  difficult  concepts  of  stability,  convergence  and  automatic 
change  of  order  and  step-size  (see  references  [2,6]).  Hachtel,  et  al. 
[49]  have  used  this  algorithm  for  electrical  network  analysis  and 
design.  Calahan  and  Orlandea  [7,8,9]  have  modified  the  Gear  algorithm 
and  used  it  in  the  ADAMS  computer  program  for  the  solution  of  the 
dynamic  system  equations.  The  Gear  algorithm  has  been  used  in  this 
dissertation  for  the  solution  of  linear  and/or  nonlinear  sets  of 
differential  and  algebraic  equations  (see  Chapter  II,  III  and  V). 

1.1.6  Matrix  Methods 

/ 

In  the  process  of  numerical  integration,  solution  of  simultaneous 
linear  equations  is  inevitable  at  every  time  step.  Theoretically,  any 
method  of  solution  by  the  L  U  factorization  procedure  (see  [6,11,44,55]) 

<v 

can  be  adopted. 

By  the  sparse  matrix  formulation  of  the  dynamical  equations, 
difficulties  withnonlinearities  in  the  differential  equations  can  be 
avoided  and  so  by  the  application  of  Gear  algorithms  and  sparse  matrix 


techniques,  large  dynamic  systems  can  be  simulated  very  effectively. 
Thus  sparse  matrix  approaches  have  been  adopted  in  this  dissertation. 

Section  2.3  of  this  dissertation  deals  briefly  with  sparse  matrix 
techniques.  Some  modifications  necessary  for  the  adjoint  analysis 
(see  Chapter  III)  are  noted  in  that  section. 

1.1.7  Thesis  Organization 

In  the  remaining  part  of  this  chapter  (i.e.,  in  Section  1.2),  a 
brief  literature  survey  on  optimization  of  dynamic  systems  is  made. 

In  Chapter  II  a  sparse  matrix  formulation  of  the  equations  of 
dynamical  systems  and  the  topics  of  stiff  integration  and  sparse  matrix 
techniques  are  discussed.  The  formulation . of  the  general  optimal  de¬ 
sign  problem  and  the  corresponding  design  sensitivity  analysis  and 
optimization  algorithms  are  described  in  Chapters  III  and  IV. 

In  the  sensitivity  analysis,  two  types  of  state  variables  occur. 

In  mechanical  problems,  the  first  type  corresponds  to  variables  like 
displacements  and  velocities  and  are  called  "Primary"  state  variables. 
The  second  type  corresponds  to  Lagrange  multipliers  of  the  constrained 
motion  and  the  spring-damper  associated  variables.  These  are  termed 
"Secondary"  state  variables  (see  Chapter  II) .  Chapter  V  deals  with  the 
organiztion  and  description  of  the  computer  program  DADS.  Numerical 
results  of  application  of  this  program  to  slider-crank  mechanism  and 
spring  reset  plow  share  mechanism  are  presented  in  Chapter  VI. 

Finally,  some  discussion,  conclusions  and  recommendations  are  given 


in  Chapter  VII. 
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1.2  Literature  Survey 

A  general  survey  of  mechanical  design  optimization  is  presented 
by  Seireg  [23].  Haug  and  Arora  [24]  have  given  a  description  of  state  4 

space  techniques  for  solving  optimal  mechanical  design  problems.  Also 

A 

a  fundamental  treatment  of  the  problem  of  optimal  design  of  constrained 
dynamic  systems  can  be  found  in  AMCP  706-192  [1].  Sevin  and  Pilkey 
[25]  have  used  the  penalty  function  technique  [26]  to  obtain  min-max 
response  of  dynamic  systems  with  incompletely  prescribed  input  func¬ 
tions.  Sevin,  Pilkey,  and  Kalinowski  [27]  have  formulated  problems  of 
optimal  design  of  mechanical  systems  subjected  to  dynamic  loads  in 
mathematical  programming  terms.  These  can  be  found  in  reference  [28]. 

They  have  studied  three  types  of  variational  problems:  (1)  Extreme 
disturbance  analysis  with  bounds  on  performance  index  for  a  given 
system  when  the  inputs  are  described  as  a  class  of  unspecified  wave¬ 
forms;  (2)  optimum  system  performance  dealing  with  bounding  a  perfor¬ 
mance  index  for  a  class  of  inputs  when  certain  system  elements  are 
unspecified  and  constraints  are  imposed  on  the  system  response;  and 
(3)  optimum  system  design  concerning  identification  of  parameters 
that  uniquely  specify  the  system,  so  that  a  performance  index  is  mini¬ 
mized  for  the  worst  disturbance  among  the  class  of  admissible  inputs. 

For  these  problems,  solution  techniques  are  based  on  linear,  nonlinear, 
and  dynamic  programming  [29]. 

Brock  [30],  Den  Hartog  [31],  Hamad  [32],  Arora,  Rim,  and  Kwak  [33], 

Afimiwala  and  Mayne  [34],  Willmart  and  Fox  [35],  McMunn  [36],  and 
Kwak  [37]  have  considered  various  aspects  of  optimization  in  vibration 


absorbers  and  vehicular  models.  Hsiao  [38]  has  considered  similar 
problems  with  "P-Norm  approximation"  and/or  "Equivalent  Functional 
Treatment"  of  cost  functional  and  performance  constraints.  Haug  and 
Arora  [39]  have  made  further  modifications  of  Hsiao's  treatments. 
However,  all  these  problems  involve  either  a  small  number  of  degrees 
of  freedom  or  they  are  relatively  easy  to  formulate  by  Newton's  laws 
(owing  to  restricted  dynamic  motions). 

In  the  field  of  intermittent  motion  of  dynamical  systems,  very 
little  work  has  been  done  so  far.  The  general  problem  of  optimization 
of  mechanical  systems  with  this  type  of  motion  has  been  treated  in 
reference  1.  Huang  [40]  and  Huang,  Haug,  and  Andrews  (41]  have 
developed  a  state-space  method  of  optimal  design  of  mechanical  systems 
with  intermittent  motion,  which  has  been  applied  to  a  cammed,  three 
mass  system.  None  of  these  methods  has  gone  into  the  consideration  of 
the  problem  of  instabilities  associated  with  a  highly  nonlinear  set  of 
differential  equations  and  with  widely  split  eigenvalues.  Thus  none 
has  taken  advantage  of  the  modern  methods  of  stiff  integration  (Gear) 
and  sparse  matrix  algorithms. 
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CHAPTER  II 

SPARSE  MATRIX  FORMULATION  OF  EQUATIONS  OF  MOTION, - 
THE  STIFF  INTEGRATION  (GEAR)  ALGORITHM, 

AND  SPARSE  MATRIX  TECHNIQUES 

2.1  Sparse  Matrix  Formulation  of  the 
Mechanical  Systems  Equations* 

2.1.1  Introduction 

The  principal  objective  of  this  section  is  to  present  a  sparsity- 
oriented  formulation  of  dynamical  equations  using  the  idea  of  writing 
the  equations  of  motion  for  each  element  (or  component)  of  a  mechanism 
(machine)  separately.  This  idea  originated  in  reference  [7]  and  is 
discussed  further  in  reference  [8].  In  this  approach,  the  Lagrangian 
formulation  of  equations  of  motion  is  adopted  and  constraint  equations 
for  various  joints  and  spring-damper  relations  are  written  separately. 
No  attempt  is  made  to  reduce  the  number  of  solution  variables  through 
the  process  of  elimination.  Thereby,  the  order  of  nonlinearity  in  the 
equations  is  kept  at  a  minimum.  The  equations  are  then  solved  numeri¬ 
cally  using  a  stiff  integration  algorithm  [2,6]  and  sparse  matrix 
[9,64,65]  techniques.  These  methods  are  discussed  briefly  in  Sections 
2.2  and  2.3. 

The  results  presented  here  are  confined  to  two-dimensional  mechani¬ 
cal  systems.  Thus,  the.  sparse  matrix  formulation  of  dynamical 

* 

Companion  reading  of  references  [7,14,63]  is  suggested  for  the 
reading  of  this  section. 
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p  3$ . 

+  l 

i=l 


(2.4a) 


and 


j  —  1,2, ...  ,k 


(2.4b) 


=  0  ,  i  =  1,2,.  ,.,p  (2.4c) 

where  are  geometrical  constraint  functions,  Q  are  generalized 
forces  (excluding  constraint  reactions),  and  p^  are  Lagrange  multi¬ 
pliers. 


2. 1.2.1  Choice  of  Coordinates 

In  two  space  dimensions,  let  0X,0Y  represent  a  set  of  coordinates 

fixed  in  an  inertial  reference  frame  and  0.x.,  O.y.  represent  a  set 

ii  ii 

of  body-fixed  rectangular  coordinate  axes  fixed  at  9^  in  the  i-th  body 

of  the  system.  Let  X^,  Y^,  <f>^  be  the  translational  and  angular  gener- 

•  •  * 

alized  coordinates  and  u^  =  X^,  v^=  Y^,  and  =  <|>^  the  corresponding 
generalized  velocities  for  the  i-th  body.  These  coordinate  systems 
are  illustrated  in  Figure  2.1.  Let 

a1  =  IVV*/ 

and 

U1  H  [Uj.Vj.wJ1 

Then,  one  can  write  the  vector  ZX(t)  of  generalized  coordinates  and 


velocities  as 
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.T  T  T 

z1  =  [u1  ,  q1  ]  =  [ui,vi,wi,Xi,Yi,4>i] 


(2-6) 


2. 1.2. 2  Spring-damper  Parameters  and  Variables 
and  Related  Equations 

To  the  system  of  equations  (2.4),  one  must  add  spring-damper 

equations.  For  the  k-th  spring-damper  pair  in  two-dimensional  systems, 

k  k 

four  variables  are  defined;  l  is  the  spring-damper  length,  v  is  the 

k  1c 

velocity  associated  with  damping  and  F  and  F  are  the  spring-damper 

A  Y 

force  components.  One  can  define  a  vector  j?k(t)  with  these  variables 
as  its  components,  i.e. 


Xk(t)  =  Uk,vk,Fk,Fk] 


(2.7) 


A  companion  vector  (t)  for  the  k-th  spring-damper  pair  con¬ 
necting  i-th  and  j-th  bodies  can  be  defined  in  the  following  way: 


(2.8) 


Figure  2.2  shows  the  spring-damper  variables  and  parameters  for  a 

spring-damper  pair.  Superscripts  identifying  the  number  of  spring- 

damper  pair  have  been  deleted  for  the  sake  of  simplicity.  The  vectors 

R.  ,R  ,r  ..,r  and  R  ..  are  position  vectors  and  s..  and  s  .  are 
i  j  sij  sji  sij  r  ij  ji 

points  of  attachment  on  the  i-th  and  j-th  bodies,  respectively.  The 
angle  a  is  measured  between  Rg„  and  the  positive  X-axis.  The  con¬ 
stants  are  spring  and  damping  coefficients,  respectively, 

k  k  k  k  k  k 

Although  K. . ,C ,v .., F  . ,  and  F  .  are  complete  notations  for 
XJ  ^-l  iiJ 

the  spring-damper  parameters  and  variables  for  the  k-th  pair 
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connecting  the  i-th  and  j-th  bodies,  superscripts,  subscripts,  or 
both  may  sometimes  be  suppressed  in  the  ensuing  discussions,  for  the 
sake  of  notational  simplicity. 

Explicit  defining  expressions  for  the  spring-damper  variables 
are  given  by  (see  [14]) 


V 


2  2 

(Uk  -  uk)  +  (Vk  -  Vk) 

J  i  J  i 


(2.9) 


(2.10) 


e  k  =  [Kk(^k  -  £k)  +  Ckvk  +  Fk]  [Uk-Uk]/.k  -  Fk  =  0 
Fx 


(2.11) 


e  k  5  [Kk(S,k  -  £k)  +  Ckvk  +  Fk]  [Vk  -  Vk]/Ak  -  Fk  =  0 

fy 

(2.12) 

lr  If  Ir  Ir 

where  (U^,V^)  and  (U^,V  )  are  the  global  coordinates  of  the 

points  of  attachment  of  the  k-th  spring-damper  pair  on  bodies  i  and 

k  k 

j,  respectively;  K  and  C  are  the  spring-damper  constants;  and 
k  k 

and  Fq  are  the  initial  length  and  constant  force  along  the  spring- 

k  k  k  k 

damper.  The  explicit  expressions  for  and  are  given  by 


k 

U.  =  X.  +  x  ..  cos  d>.  -  y  ..  sin  * 

l  l  sij  i  sij  i 


(2.13) 


k 

V .  =  Y .  +  x  .  .  sin  <b .  +  v  .  .  cos  d> . 

1  1  sij  1  '  sij  1 


k 

U .  =  X .  +  x  .  .  cos  <j> .  -  y  .  .  sin  A . 

J  J  sji  j  sji  j 


(2.14) 


(2.15) 
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Vk  =  Y.  +  x  .  .  sin  tj).  +  y  .  .  cos  d>.  (2.16) 

J  J  sji  1  sji 


where  (x  . . , 
si  3 

s . .  and  s  .  . , 
H  Ji 

axes. 


y  ..)  and  (x  ..,y  ..)  are  the  coordinates  of  the  points 
sij  sji 

with  respect  to  their  respective  body-fixed  coordinate 


The  functional  symbols  e  ,  e  ,  ,e  ,  ,e  ,  in  Eqs.  (2. 9) -(2. 12)  are 

F  Vk  Fk  F* 

written  to  indicate  the  small  values  the  expressions  will  take  during 
Newton  iterations  in  the  dynamic  analysis  (see  Section  2.2  and 
reference  [6]). 


2. 1.2. 3  Constraint  Equations 

In  two  dimensions  there  are  two  principal  types  of  joints; 

revolute  and  translational.  Figure  2.3  shows  parameters  and 

coordinates  associated  with  a  revolute  joint.  In  the  figure,  r^  is 

the  position  vector  of  a  point  P_.^  on  body  j  with  respect  to  a  point 

P. ..on  body  i.  When  r  =  0,  P..  and  P..  coincide  and  they  define  a  re- 
ij  P  ij 

volute  joint.  The  loop  closure  relation  of  the  position  vectors 
gives . 


R.  +  r .  .  +  r  -r.  .  -  R.  =  0 
1  13  P  3 1  3 


(2.17) 


Let  (x.,y.)  and  (x.,y.)'be  the  coordinates  of  P..  and  P..  respectively 
i  i  3  3  13  Ji 

with  respect  to  body-fixed  axes.  With  r^  =  0,  the  constraint  equa¬ 
tions  for  a  revolute  joint  between  bodies  i  and  j  are  obtained  from 


Eq.  (2.17)  as, 
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<i>  =  X.  +  x.  cos  4> .  —  y  •  sin  <j> .  -  X.  -  x  .  cos  d> .  +  y .  sin  <f> .  =  0 

X  i  x  1  i  i  j  j  Tj  Tj 


(2.18) 

and 


$>  =  Y .  +  x .  sin  <j> .  +  y .  cos  6.  -  Y.  -  x.  sin  d> .  -  y .  cos  d> .  =  0 

Y  i  i  Ti  i  Ti  J  1  j  1  TJ 

(2.19) 

Figure  2.4  shows  parameters  and  coordinates  associated  with  a 

i  i 

translational  joint.  Here  P..  and  P..  are  the  points  of  intersection 

3 1 

of  perpendiculars  drawn  from  the  origins  (centers  of  mass)  0^  and  0_. 

onto  a  straight  line  that  is  parallel  to  the  line  of  relative  motion 

between  the  bodies.  The  vectors  r..,r..,r  have  the  same  meaning  as 

ij  Ji  P 

before  and  6.  and  6.  are  angles  between  the  vectors  r;.  and  r..  with 
13  13  Ji 

the  body-fixed  axes  O.jX^  and  0 ,  respectively.  The  loop  closure 
condition  of  the  position  vectors  again  gives,  after  elimination  of 
r^,  the  following  constraint  equations  for  a  translational  joint 
between  bodies  i  and  j  (see  [14,63]) 

n  x  i  x  i  ri  i 

-  X.  cos(<j>.  +  6.)  -  Y.  sin(<j>.  +  6.)  -  Jx4.  +  y^  =0 

J  J  J  3  33  33 

(2.20) 

*  =  <fr.  +  6.  -  -  6.  -  0  (2.21) 

<f>  i  i  3  3 

In  Eqs.  (2.18)  to  (2.21),  the  coordinates  (x^,y^)  on  body  i  and 
(x ^  » y j )  on  body  j  (hence  the  derived  parameters  6^  and  6.)  depend  on 
design  parameters  that  define  the  geometry  of  the  bodies. 


/2  2 

'x .  +  y . 
i  i 
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2. 1.2. 4  Primary  and  Secondary  State  Variables, 

Design  Parameters 

In  the  problem  of  optimal  design  of  constrained  dynamic  systems, 
state  variables  and  design  parameters  are  encountered  [1].  In  the 
present  formulation,  two  kinds  of  state  variables  are  defined: 

Primary  State  Variables:  Variables  whose  time  derivatives  appear 
in  the  equtions  of  motion  or  in  related  equations  are  called 
"Primary  State  Variables". 


iT  iT 

For  the  planar  systems  treated  here,  u  ,q  ,  and  1^,1 


J1  J2 


A 


(the  last  k  terms  being  the  lengths  of  the  k  spring-damper  pairs 


connecting  i-th  body)  are  Primary  State  Variables  related  to  the 
* 

i-th  body. 

Secondary  State  Variables:  Variables  that  appear  without  their 
time  derivatives  (i.e.,  appear  algebraically  only) in  the  equations  of 
motion  and  related  equations  are  called  "Secondary  State  Variables". 

Let  h1  =  h1 (t)  denote  the  Lagrange  multiplier  vector  correspond¬ 
ing  to  the  constraints  on  the  i-th  body.  Then  the  variables 


J  v  do  ^k 

v  ,F  ,  F  ,...,F  ,F  ,F  ,...FY  are  secondary  state 

A  A  X.  Y  X 


H  ,v  ,v  ,  . 
variables  related  to  the  i-th  body. 

The  vector  b  denotes  the  design  parameter  vector  for  the  entire 
system. 


2. 1.2, 5  Element  State  Equations  of  Motion 

In  the  Lagrangian  formulation  [66,67,68]  of  the  equations  of 
motion  of  a  two  dimensional  constrained  mechanical  system,  when  all 

*  k  .  k 

l  is  kept  with  other  spring-damper  variables  in  X  (t)  for  advantages 

in  computer  programming. 
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the  generalized  coordinates  and  Lagrange  multipliers  are  taken  to 
be  independent,  the  state  equations  of  motion  for  the  i-th  body  of  the 
system  can  be  written  as  (cf.  Eqs.  (2.1)  and  (2.4)), 

Pi(b)z1  +  f  =  0  (2.22) 

where 

pV)  = 

M.  and  J.  are  the  mass  and  moment  of  inertia  of  the  i-th  body, 
li 


M. 


M. 


J. 

l 


-1 


-1 


-1 


(2.23) 


y,,k  =  l,2,...,p.,  are  the  components  of  the  Lagrange  multiplier  vector 

K.  1 

yX(t),  and  Q„  ,  Q__  ,Q,  are  generalized  forces  (including  the  contri- 
—  X .  Y .  <p . 

1  1  yl 

butions  from  spring-dampers).  Equation  (2.23)  indicates  that  and 

JL  may  be  taken  as  design  parameters.  Note  that  the  last  three  equa- 

•  i 

tions  of  Eq.  (2.22)  contain  3  explicitly.  This  structure  is  intro¬ 
duced,  intentionally,  in  order  to  increase  sparsity  (see  Section  2.3). 


2. 1.2.6  Element  State  Equations  for  a 
Slider  Crank  Mechanism 

To  illustrate  the  procedure  for  writing  the  system  equations  of 
motion,  a  slider  crank  mechanism  is  considered  in  this  subsection. 

The  radial  slider-crank  mechanism  is  a  complex  of  rigid  bodies  that 
move  in  a  plane.  Figure  2.5  shows  the  approximate  initial  position  of 
such  a  mechanism.  Link  1  is  ground ,  link  2  is  the  crank  shaft,  link  3 
is  the  connecting  rod  or  coupler,  and  link  4  is  the  piston  (or  slider). 
A  spring  damper  pair  is  attached  between  link  4  and  ground  (Figure  2.5). 
There  is  a  translational  joint  (type  2)  between  bodies  4  and  1  and  a 
revolute  joint  (type  1)  between  each  of  the  following  pairs  of  bodies: 
1  and  2,  2  and  3,  and  3  and  4.  The  figure  also  indicates  the  follow¬ 
ing  design  parameters: 

b^  =  The  spring  constant  of  the  spring 

b^  =  Height  of  the  points  of  attachment  of  the  spring 

b^  =  Half  of  the  length  of  the  uniform  coupler. 

Gravitational  forces  are  excluded  from  the  present  simulation  of  the 
mechanism.  To  illustrate  the  use  of  the  present  technique,  only  the 
equations  of  motion  and  constraint  involving  body  4  are  written. 

For  body  4,  there  will  be  4  constraint  equations  corresponding 
to  (revolute)  joint  3  and  (translational)  joint  4.  They  can  be 
written  as  (using  notations  of  previous  subsections). 


<S>  =  X„  +  x„  cos  4>_  -  v„  sin  <)>„  -  X.  -  x,  cos  ct>,  +  y,  sin  <f>  =0 

A.  _>  j  J'J  J  4  4  44  4 


(2.25  cont.) 


i-aapTis 
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(xg^,y  ^),  (xs4i>ys4]_)  are  respectively  the  coordinates  of  the  points 

and  of  attachment  of  the  spring-damper  pair. 

In  this  problem,  contributions  to  the  generalized  forces 

Qv  »  Qv  >  end  Q  come  only  from  the  spring-damper  forces.  Thus, 

X4  Y4  ^4 


x  -  v-  -v 

4  1 


Q  =  F  (=  -Q  ) 
*4  1  1 


(2.28) 


\  =  “  FX(xs41  Sin  ^4  +  ys41  C0S  V  +  FY(xs41cosn-ys41Sin  V 


From  Eqs.  (2.22)  to  (2.24),  the  state  equations  of  motion  for  the 
4-th  body  can  be  written  as. 


• 

r~  v 

'  — — 

M4 

U4 

\  3X4  "k 

0 

• 

t  3ik 

M4  0 

v4 

QY.  , X  dY.  pk 

0 

4  k=l  4 

• 

r  3*k 

J4 

w4 

^4  k=l  3 ^4  Pk 

0 

-1 

*4 

‘  u4 

0 

(o 

i 

M 

*4 

"  V4 

0 

-1 

*4 

'  W4 

0 

(2.29) 
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The  equations  of  motion  and  constraint  relations  for  other  bodies 
of  the  mechanism  can  be  written  in  a  similar  manner.  However,  manual 
formulation  of  the  equations  is  unnecessary,  since  the  computer  program 
ADAMS-2D  [14]  automatically  generates  the  necessary  code  for  all  the 
expressions  and  equations.  Only  data  defining  parameters  such  as 
mass,  moments  of  inertia,  locations  of  centers  of  mass,  location  of 
joints  and  points  of  attachment  of  spring-damper  pairs,  and  spring- 
damper  constants,  (see  [14,63])  are  to  be  provided  by  the  user. 


The  global  vector  of  Lagrange  multipliers  may  be  denoted  by  p  which, 
with  its  components  H^jI^*****  is 


y  =  [y1,y2> 


2m-l,M2mJ 


(2.31) 


The  global  vector  of  spring-damper  variables  may  be  denoted  by 
which  with  its  components  1^,1^, . . . , . is  given  by 


"  T  T  T~l^ 

1  2  s 

l  =  t  ,L 


, &2  >  •  • 


,1. 

4s-l  4s 


(2.32) 
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Now  the  state  equations  of  motion  for  the  entire  system  can  be 
written  in  the  form  (deleting  for  simplicity  the  underlines  of  the 
vector  variables) 


F(t,z,z, £,p,b)  =  P(b)z  +  f(t,z,£,p,  b)  =  0 


(2.33) 


where 


P(b) 


PX(b)  0 

P2(b)_ 

0  ‘  pn(n) 

(6n  x  6n) 


(2.34) 


f 


T  T 
f  f 
1’  2’ 


From  equations  (2.18)  to  (2.21)  one  may  summarize  the  constraint 
equations  in  the  form 


$(z,b)  =  0 


(2.35) 


From  equations  (2.9)  to  (2.12),  one  may  write  the  equations 
related  to  spring-damper  pairs  as 


\  =  IS.  +  £(z,£,b)  =  0 


(2.36) 


where 


£  is  the  algebraic  part  of  5,  and  I  is  defined  as 


(2.37) 
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I  = 


0  0 
1  0 


0  0 
1  0 


(2.38) 


(4s  x  4s) 


The  system  of  equations  (2.33),  (2.35),  and  (2.36)  is  solved  by 
Gear's  predictor-corrector  algorithm  and  sparse  matrix  techniques 
discussed  briefly  in  Sections  2.2  and  2.3. 

2.2  Stiff  Integration  (Gear)  Algorithm 

2.2.1  Introduction 

Any  system  of  differential  equations  that  has  widely  split 
eigenvalues,  at  least  locally,  is  called  a  stiff  system.  If  a 
nonlinear  system  N  is  approximated  by  a  linear  system  L  around  some 
point  t,  the  eigenvalues  of  L  are  the  local  eigenvalues  of  N  at  t 
with  respect  to  L.  The  mechanical  system  of  Fig.  2.6  represents  a 
stiff  mechanical  system.  As  indicated  in  Chapter  I,  mechanical 
systems  like  the  ones  discussed  in  Chapter  VI  may  not  be  stiff  ini¬ 
tially,  but  they  may  unpredictably  become  stiff.  It  has  been  shown 
in  references  [2,6]  that  neither  Runge-Kutta  nor  Adams-Bashforth 
nor  Adams-Moulton  algorithms  are  suitable  for  the  solutions  of  such 
systems,  for  stability  reasons  [2,6]. 


=1000 


Figure  2.6  Stiff  Mechanical  System. 
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In  1969,  Gear  developed  stiffly  stable  [2,6,47,48]  multistep 
algorithms,  which  are  well  suited  for  the  solution  of  stiff  systems. 
Originally  he  considered  systems  of  the  form  of  Eq.  (1.3)  and  employed 
his  criteria  of  stiff  stability  to  derive  the  algorithm.  Later  he 
showed  [46]  that  the  same  algorithm  can  be  used  for  mixed  systems  of 
differential  and  algebraic  equations.  One  must  understand  the  techni¬ 
cal  details  of  the  subject  of  multistep  numerical  predictor-corrector 
algorithms  and  the  concepts  of  stability  and  convergence  and  automatic 
change  of  order  and  step-size  in  order  to  have  proper  command  of  the 
Gear  algorithms  and  their  implementation  in  computer  program  DIFSUB 
(see  [7]).  These  are,  however,  availabe  in  standard  references  [2,6]. 
and  will  not  be  treated  here.  The  basic  idea  of  automatic  control  of 
order  and  step  size  can  be  stated  briefly  as  follows: 

Let  t  =  t  represent  the  current  time  instant  and  h  and  k  be  the 
n 

current  step  size  and  order  of  the  numerical  integration.  Let  be 

the  local  truncation  error  [6]  and  e  the  maximum  allowable  e_. 

max  1 

The  basic  step  control  algorithm  is  then  to  execute  one  step  and 
test  whether  the  relation 


£ 


T 


< 


e 

max 


(2.39) 


is  satisfied.  If  it  is,  the  step  is  accepted.  Otherwise,  the  step 
is  rejected  and  a  smaller  step  size  h  =  ah(a  <  1)  is  used.  The 
exact  step  size  to  use  for  executing  the  next  step,  or  for  repeating 
the  rejected  step,  is  given  by  the  choice  of  a  as  the  maximum  value 
computed  from  three  expressions  [6]  for  local  truncation  errors 
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corresponding  to  the  orders  k,  k  -  1,  and  k  +  1.  The  maximum  of  the- 
three  ex's  gives  the  maximum  allowable  step  size  and  the  corresponding 
order  is  the  optimum  order  to  compute  values  at  t  =  t 

2.2.2  The  Gear  Algorithm  and  the  Mixed  System 
of  Differential  and  Algebraic  Equations 

The  k-th  order  Gear  algorithm  for  the  solution  of  the  mixed 
system  of  differential  and  algebraic  equations  (1.1)  and  (1.2)  is  an 
implicit  algorithm  of  the  form 


-n+1 


V  -n+1  -  (Vn  +  al?n-l  + 


ak-lYn-k+l^ 


k-1 

1  a.  y 

j=o  J  n_:J 


(2.40) 


where  h  is  the  time  step,  y  ,y  , ...  are  the  values  of  y  at  time 

instants  t  . -,t  and  a„»a, , .  •  •  »a,  ,  ,B„  are  the  (k  +  1)  coeffi- 

n+1  n  0  1  k-1  0 

cients  known  as  Gear  coefficients  for  this  multistep  algorithm  (2,6]. 

One  proceeds  from  the  n^  to  the  (n  +  l)St  time  step  by  solving 
Eq.  (2.40),  together  with  Eq.  (1.1)  at  t  =  t  i.e. 


-  o 


(2.41) 


Linearization  of  Eq.  (2.41)  gives  the  Newton  formula  [6]  at  the 
s  t 

(n  +  1)  time  step: 


3F(n) 

3Z 


SF(m) 
+  - 

3y 


(2.42) 
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where 

.  (m)  (nrKL)  (m) 
A^  =y  -  y 

•  (m)  •  (nH-1)  •  (m) 

=  y  -  y 


(2.43) 


m  being  the  iteration  number  in  the  Newton  method  of  solving  the 
algebraic  equations.  The  time  step  counter  n  has  been  dropped,  for 
simplicity  of  notation. 


k-1 


Since  7  a.y  .  remains  invariant  for  Newton's  iteration  at 
j=0  2~n~2 


st 


the  (n  +  1)  time  step,  one  obtains  from  Eq.  (2.40) 


a/")  -  -eoh  Ay<m) 


or 


AZ 


(m) 


1  .  (m) 

VA2 


(2.44) 


Hence,  Eq.(2.42)  becomes  (after  the  substitution  of  the  second  of 
Eqs.  (2.44)) 


3F<"> 

By 


9F(m) 

3Z 


(2.45) 


s  t 

This  is  called  the  "corrector  formula"  for  y  at  the  (n  +  1)  time 
step. 


Equation  (2.45)  together  with  the  second  equation  of  Eq.  (2.44), 
updates  both  y  and  y,  which  are  required  in  The  iteration  is 

continued  until  the  right-hand  side  of  Eq.  (2.45)  is  less  than  a  pre¬ 
assigned  small  quantity.  Since  the  Jacobian  matrix  on  the  left-hand 
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side  of  Eq.  (2.45)  is  of  the  same  structure  for  each  iteration,  this 
procedure  matches  ideally  the  requirements  for  code  generation  for 
sparse  matrix  algorithms  discussed  in  Section  2.3. 

It  should  be  noted  that  in  this  procedure  even  a  linear  differen¬ 
tial  equation  will  generally  require  more  than  one  Newton  iteration  for 
corrector  convergence.  Moreover,  this  procedure  may  be  adopted  for 
any  implicit  algorithm. 

The  above  procedure  for  the  evaluation  of  y  , ,  and  y  ,  is  due 

-n+1  -n+1 

to  Calahan  and  Orlandea  [7,8,9]  and  has  been  used  in  the  ADAMS 
program  [ 7 ] . 

For  various  advantages  in  computer  programming,  implicit  multi- 

step  predictor-corrector  algorithms  are  recast  into  canonical  matrix 

representations  (see  references  [2,6,7]).  For  that  purpose,  the  Gear 

algorithm  is  recast  into  such  a  representation  with  the  help  of 

Nordsieck  vector  z  [2,6,7]  that  is  defined  as 
- -n 

'  2  "  k  M  T 

z  5  [yT,’hyn»h  y  /2!,...,h  y„  /k!]  (2:46) 

— n  n  n  n  n 


where 


*  »  v  " 
yn  ,yn  ’ 


are  the  1st,  2nd,  ...,  kC^  derivative  respectively  of  a  single  compo¬ 
nent  y  at  t  =  t  .  The  Nordsieck  array  z  is  defined  bv 
n  - ^  ~n 


z 

~n 


[ y  ♦  hyn ' 
-n  -  n 


>h 


V/2! 


,hkZ<k>/k!] 


r 


* 


(2.47) 


All  the  updated  values  of  the  Nordsieck  vectors  are  required  for  pre¬ 
diction  of  the  solution  variables  that  are  used  as  initial  estimates  in 
the  corrector  iterations.  It  should  be  noted  here  that  Eqs.  (2.44)  and 
(2.45)  give  only  the  first  two  components  of  the  Nordsieck  vectors. 

Orlandea  [7]  has  shown  that  all  the  components  of  the  Nordsieck 
vectors  can  be  obtained  from  the  corrector  iteration  formulas  of  the 
form: 

3F(m)  1  3F(m) 

3y  6Qh  3y 

where  the  time  step  counter  n  has  again  been  suppressed  for  simplicity, 

the  vector  represents  the  vector  of  the  i-th  components  of 

at  the  m-th  Newton  iteration,  and  c  .,  i  =  l,2,...,k+l,  are  the 

zx 

coefficients  of  the  transformed  Gear  algorithm.  Their  values  are 
given  in  Table  2.1  (see  reference  [6]). 

The  present  form  of  DIFSUB,  however,  does  not  iterate  for  the 
values  of  the  Nordsieck  vector  components  other  than  the  first  two. 

They  are  evaluated  from  Eq.  (2.48)  after  corrector  convergence  for 
the  first  two  components. 

2.2.3  Starting  of  Multistep  Algorithms 
Multistep  numerical  algorithms  are  not  self  starting.  Generally 
a  single  step  algorithm  is  used  at  least  k  times  before  a  multistep 
algorithm  can  be  initiated. 

In  the  ADAMS  program.  The  Gear  algorithm  is  implemented  through 
the  subroutine  DIFSUB  and  initially  the  order  is  taken  to  be  1.  The 


Az 


(m)(i)  _  _  Czi  F(m) 
Czl  " 


(2.48) 


Table  2.1 


Coefficients  of  Stiffly  Stable  Methods 
in  Canonical  For:?. 


k 

2 

3 

4 

5 

6 

2 

6 

25 

120 

720 

c  , 
zl 

3 

11 

50 

274 

1764 

3 

11 

50 

274 

1764 

Cz2 

3 

11 

50 

274 

1764 

1 

6 

35 

225 

1624 

Cz3 

3 

11 

50 

274 

1764 

1 

10 

85 

735 

Cz4 

11 

50 

274 

1764 

1 

15 

175 

Cz5 

50 

274 

1764 

1 

21 

cz6 

274 

1764 

c  i 

z7 

1 

1764 

first  order  Gear  algorithm,  being  exactly  the  backward  Euler  algorithm, 
is  itself  a  single  step  algorithm  that  serves  the  purpose  of  initial¬ 
ization. 


2.2.4  Corrector  Formulas  for  the  Dynamical 
System  Equations 

For  the  state  equations  of  motion  gi,ren  by  Eqs.  (2.33),  together 
with  the  constraint  equations  (2.35)  and  the  spring-damper  relations 


(2.36),  the  corrector  formula  (2.45)  can  be  written  as  (deleting 
underlines  and  superscripts) , 


(2.49) 


where  the  following  relations  have  been  used: 


3F  =  3f  ll  «  li 

3z  3z  ’  3z  3z 


3F  =  If  ll  =  li 

3y  ~  3w  ’  32,  35, 


(2.50) 


3F  =  3f 
32,  =  32, 


i  T 

Moreover,  when  z  =  [X^ Y±,  ^’“i^i^i5  instead  of 
[ui,Vi’Wi,Xi’Yi^i]T»  one  °^ta:i-ns  (see  Eq.  (2.24)), 


(2.51) 


In  that  special  case,  Eq.  (2.49)  can  be  written  as 
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(2.52) 

A  descriptive  tableau  form  for  the  corrector  formula  of  Eq.  (2.49)  can 
be  given  as 


Body  1 


Body  n 


constraint 
related 
sparse 
submatrices 
[with  elements 

lT 

!  of 


constraint  related 
sparse  submatrices  with 
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elements  of  — 

3q 
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hamper 
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submatrices 
with  elements 
3F 

°f  a 


spring-damper  related  submatrices  of  the  form: 
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3u 


ii 

3q 


ii 
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V 

ta2j 

Ay 

-5 

AZ 

(2.53) 


The  matrix  in  the  left-hand  side  of  this  corrector  formula  is 
known  as  the  corrector  Jacobian  matrix.  A  detailed  description  of 
the  nonzero  positions  in  the  Jacobian  matrix  of  such  a  tableau  for 


the  four-bar  slider-crank  mechanism  of  Fig.  2.5  is  given  by  Fig.  2.7. 
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Figure  2.7  Symbolic  Listing  of  the  Nonzero 
Entries  in  the  Jacobian  Matrix 
for  the  Example  Slider-Crank 
Mechanism. 
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2.3  Sparse  Matrix  Techniques 

2.3.1  Introduction 

In  order  to  solve  the  system  of  simultaneous  linear  equations  of 
Eq.  (2.52),  some  matrix  method  must  be  employed.  Sparse  matrix  tech¬ 
niques  enhance  speed  of  computation  in  such  systems.  These  techniques 
are  described  briefly  in  the  following.  If  less  than  30%  of  the 

entries  in  a  square  matrix  A  are  nonzero,  the  matrix  A  is  considered 

~nxn 

to  be  sparse  and  storing  the  matrix  as  a  two-dimensional  array  becomes 
inefficient.  Consideration  of  matrix  sparsity  is  extremely  important 
for  speed  of  computation  in  problems  of  mechanical  system  analysis 
(particularly  dynamic  systems  considered  later  in  this  dissertation) . 
This  consideration  outweighs  the  difficulties  encountered  in  solving 
a  large  set  of  simultaneous  linear  equations  to  which  many  physical 
systems  can  be  ultimately  reduced  [50,51]. 

Sparse  systems  of  simultaneous  linear  equations  generally  result 
from  the  solution  of  the  following  classes  of  equations: 

(1)  Ordinary  differential  equations  and/or  algebraic  equations,  where 
after  time  discretization  and/or  linearization  an  irregularly 
structured  matrix  is  encountered  (cf.  Fig.  2.6); 

(2)  Partial  differential  equations,  where  after  discretization  of  the 
spatial  variables  by  finite  difference  or  finite  element  tech¬ 
niques  [52,53,54]  a  regularly  structured  matrix  is  handled. 

This  research  is  concerned  only  with  problems  of  the  first  class 
(as  is  evident  from  the  last  section  and  the  developments  in  the 


subsequent  chapters) . 
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In  sparse  matrix  methods,  operations  involving  zeros  are  avoided 
by  using  structural  information  (the  positions  of  nonzero  entries) 
that  is  stored  in  a  compacted  form.  One  way  of  doing  this  is  to  store 
the  row  and  column  indices  of  each  nonzero  element  in  two  vectors  I 
and  J  and  the  value  of  the  elements  in  a  third  vector  G.  This  method 
is  called  "i-j"  ordering.  According  to  Calahan  [10]  this  is  the  most 
convenient  method  and  can  be  easily  converted  to  other  methods  of 
compacting  the  data, such  as  the  threaded  list  method  and  the  bit  map 
method,  which  are  discussed  in  detail  in  [10]. 

2.3.2  Solution  of  Simultaneous  Linear 
Algebraic  Equations 

There  are  two  general  methods  for  solving  a  set  of  simultaneous 
linear  algebraic  equations 

Ax  =  B  (2.54) 

These  are  (1)  the  Gaussian  elimination  method  and 

(2)  L  U  factorization  method  [6,11,44].  Although  theoretically  they 
are  somewhat  interconnected,  the  solution  techniques  involve  two 
different  algorithms.  The  L  U  factorization  method  is  preferable  for 
sparse  matrix  techniques.  Since  the  inverse  of  a  sparse  matrix  may 
be  full,  whereas  L  U  factors  may  retain  sparsity  [10],  the  number  of 
operations  in  the  first  method  is  much  larger  than  that  in  the  second 
[6].  The  method  of  L  U  factorization  is  now  briefly  described. 
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Let  L  and  U  be  of  the  form: 


Then  to  get 


hi’ "hi 


hi’ "hi' 


1  •  • *  u~  .  • •  •  u 
•  lj  In 


(2.55) 


A  =  L  U 


(2.56) 


one  has  with  a„  as  the  elements  of  A, 


4-n  “  a-!-!  "  l  ^ik  ukj  ’ 


ij  ij 


i  >  j 


(2.57) 


u«T«“k A1  ’ 


i  <  j 


(2.58) 


which  determine  the  matrices  L  and  U  for  the  matrix  4(nxn)- 

This  method  is  employed  in  the  Crout  algorithm  [55]  for  L,  U  factor¬ 


ization. 


After  L  U  factorization  of  A,  the  solution  of  Eq.  (2.54)  is 
obtained  as  follows: 


L  U  x  =  B 


(2.59) 
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That  is,  the  number  of  fills  should  be  kept  to  a  minimum.  This  can  be 
done  by  a  suitable  permutation  of  rows  and  columns.  Such  an  operation 
is  known  as  Optimal  Ordering  (or  pivoting) .  Some  simple  examples 
demonstrating  the  effectiveness  of  such  ordering  can  be  found  in  [6,7]. 
Over  and  above  the  solution  efficiency  and  the  minimization  of  the 
size  of  the  generated  code,  there  is  another  important  reason  for 
optimal  ordering.  Most  of  the  equations  in  physical  problems  are 
nonlinear  and  the  entries  in  the  matrix  vary  from  step  to  step  of  a 
solution  process.  In  such  cases  optimal  ordering  prevents  generation 
of  zero-valued  pivots  and  costly  regeneration  of  the  solution  code. 

The  detailed  discussions  of  the  subject  of  optimal  ordering  is 
beyond  the  scope  of  this  dissertation.  They  can  be  found  in  refer¬ 
ences  [56,57,58,59].  Codes  such  as  OPTORD  [60,49]  and  MOOP  [61]  are 
two  of  several  computer  programs  that  can  be  used  for  optimal  ordering. 
Although  optimal  ordering  aims  at  the  largest  (absolute  value)  non¬ 
zero,  row-column  entry  as  the  pivot  and  minimization  of  the  number  of 
fills,  generally  both  cannot  be  achieved  simultaneously  and  the 
algorithm  chooses  from  among  the  larger  than  average  pivots,  the 
one  which  results  in  the  minimum  number  of  fills.  The  mode  of 
operation  for  optimal  ordering  can  be  visualized  from  the  following 
considerations . 

Let  K  be  an  integer  between  1  and  n  that  represents  the 
step  associated  with  the  permutation  and  L  U  factorization  of  the 
remaining  matrix  at  that  pivot  step.  Every  nonzero  element  in  the 
residual  submatrix  of  dimension  (rt  -  K  +  1)  is  considered  as  a 
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candidate  for  the  pivotal  element.  Associated  with  each  nonzero 
element  A(I,J)  in  the  submatrix  is  a  weighting  function: 

W(I,J)  =  (LROWI  -  1)  (LCOLJ  -  1) 

where  LROWI  is  the  number  of  nonzero  elements  in  the  row  containing 
A(I,J)  and  LCOLJ  is  the  number  of  nonzero  elements  in  the  column  con¬ 
taining  A(I,J).  The  term  W(I,J)  represents  the  number  of  multiplica¬ 
tions  required  if  A(I,J)  were  selected  as  pivotal  element.  For  each 
row  with  I  fixed,  W(I,J)  is  determined  only  for  elements  that  are 
numerically  acceptable;  i.e.,  if 

|A(I,J)|  >  e  •  AVG(I) 

where  0  <  e  <  1  is  a  user-specified  tolerance, 

i  N* 

AVG(I)  =  ~  l  |A(I,J)  | 

I  J=1 

and  N^.  is  the  number  of  nonzero  elements  in  row  I.  The  algorithm 
attempts,  in  a  systematic  manner,  to  locate  the  element  with  the 
largest  absolute  value  and  the  smallest  value  for  the  weighting 
function. 


2.3.4  Column  Ordering  of  a  Matrix 
It  has  been  mentioned  in  Subsection  2.3.1  that  the  nonzero 
values  of  the  original  matrix  are  stored  in  a  vector  G.  They  can  be 
stored  either  in  row-wise  order,  taking  one  row  after  another,  or  in 
column-wise  order,  taking  one  column  after  another.  The  optimal 
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ordering  algorithms  can  handle  both  representations.  However,  the 
subroutines  VMSP,  VMNP,  and  VMBP,  used  in  ADAMS  programs  for  L  U 
factorization  of  matrices  and  numerical  solution  of  linear  equations, 
consider  only  the  column  order  representation.  The  various  aspects 
of  all  these  subroutines  are  discussed  in  reference  [62],  Many  other 
subroutines  performing  similar  functions  can  be  found  in  [61]. 

In  the  original  arbitrary  "i-j"  ordering,  the  row  and  column 

indices  are  stored  in  vectors  NPOSR  and  NPOSC,  respectively,  and  the 

nonzero  elements  of  the  matrix  A  are  stored  in  the  vector  G.  To 

transform  the  "i-j"  ordering  into  column  ordering,  a  permutation 

vector  is  generated  to  cause  the  elements  a.,  of  A  to  be  entered  into 

ij 

the  appropriate  locations  in  the  vector  G. 

The  process  of  converting  the  code  from  "i-j"  ordering  to  column 
ordering  is  illustrated  in  Table  2.2.  A  unique  number  NRANK(I)  is 
assigned  to  each  of  the  NG  nonzero  entries  in  the  matrix,  in  such  a 
way  that  NRANK  increases  as  one  proceeds  down  a  column.  No  value  of 
NRANK  for  a  certain  column  is  greater  than  that  for  any  other  column 
to  its  right.  After  the  initialization  of  a  counter  NCOUNT  in  Step  2 
and  reordering  the  NRANK  and  NCOUNT  in  Step  3,  Step  4  generates  the 
row  indices  IA(I)  of  the  desired  column  ordered  matrix  and  Step  5 
generates  the  points  JA(I)  to  the  subscripts  of  the  indices  of  the 
first  element  in  each  column.  Step  6  generates  the  permutation  vector 
NPOS  for  column  ordering  of  the  original  G  vector.  This  algorithm  is 
implemented  in  subroutine  S08000  of  ADAMS  2-D  [14]  and  SEP  of  ADAMS 
3-D  [7].  Simple  numerical  examples  of  this  procedure  can  be  found  in 
reference  [63]. 


Table  2.2 


Step  1 

Step  2 
Step  3 

Step  4 


Step  5 

Step  6 


Column  Ordering  of  a  Random  Matrix 


N  =  dimensions  of  matrix 
Calculate  the  rank  of  each  nonzero  entry: 
NRANK(I)  =  NPOSR(I)  +  (N  + 1) *NPOSC (I) ,  1=1, NG 
Initialize  a  vector  NCOUNT(I)  =  I,  1=1, NG 
Reorder  NRANK  from  the  smallest  to  the  largest. 
Reorder  NCOUNT  along  with  NRANK 
Generate  a  vector  IA(I),  1=1, NG,  that  gives  the 
row  indices  of  the  permuted »G  vector: 

IA(I)  =  NRANK(I)  -  (NRANK(I)/ (N+1))*(N+1) 

Generate  a  vector  JA(J) ,  J=1,N,  that  points  to 
the  index  I  of  IA(I)  of  the  first  nonzero  entry 
in  each  column,  and  JA(N+1)  =  NG  +  1. 

Generate  a  vector  NPOS (NCOUNT (I) )  =  I,  1=1, NG. 


52 


In  practice,  one  makes  the  assignments  G(NPOS(I))  =  a^_.  to  put 
G(I)  in  column  order.  The  vector  NPOS  is  generated  only  once  and  is 
never  modified  during  the  execution  of  the  program.  With  the  help  of 
this  vector,  any  element  of  G  can  be  directly  accessed  and  modified, 
if  desired. 

In  the  sensitivity  analysis  (see  Chapter  III  of  this  report), 
the  solution  of  adjoint  equations  requires  column  order  representation 
of  the  transpose  of  the  original  matrix.  To  achieve  this  the  original 
vectors  NPOSR  and  NPOSC  are  stored  in  another  set  of  vectors  NP0SC1 
and  NP0SR1,  respectively,  and  subroutine  S08000  is  called  to  generate 
the  vectors  NPOS  and  JP1  corresponding  to  original  vectors  NPOS  and  JA. 
This  is  done  in  subroutine  DYNANL  (see  Chapter  V) . 


2.3.5  Matrix  Vectorization 

Consider  a  column  of  a  sparse  matrix  having  the  nonzero  row 


positions  shown  in  Fig.  2.8.  This  structure  is  described  in  the 


r1-*1  column 


Figure  2.8  Matrix  Vectorization. 
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conventional  ordered  list  as 

25,  26,  27,  28,  32,  35,  36,  37,  38  (2.64) 

A  list  enumerating  all  the  row  positions  in  the  column  is  called 
"scalar  storage".  For  large  matrices,  additional  savings  in  computer 
memory  and  execution  time  can  be  realized  by  further  compacting  the 
column-ordered  code.  This  is  done  by  a  subroutine  VECTOR  after  the 
optimal  ordering  operation.  The  algorithm  obeys  the  following  rules. 

If  in  a  given  column,  there  is  a  set  of  two  or  more  contiguous 
row  positions,  only  the  first  and  last  row  indices  are  retained  in  the 
vector  IA  of  Table  2.2.  This  implies  that  all  elements  in  between 
them,  the  boundary  terms  inclusive,  are  nonzero. 

On  the  other  hand,  if  there  is  an  isolated  element  in  a  column, 
with  no  nonzero  adjacent  term,  its  row  index  is  set  negative. 

JA  of  Table  2.2  is  updated  to  reflect  the  changes  in  IA  and  has 
the  same  meaning. 

til 

Thus  after  vectorization,  the  r  column  shown  in  Fig.  2.8  is 
described  in  the  ordered  list  as 

25,  28,  -32,  35,  38  (2.65) 

Given  the  vectorized  and  permuted  description  of  the  nonsingular 
matrix,  subroutine  VMSP  (Vectorized  Matrix  Symbolic  Preprocessor)  can 
be  used  to  generate  a  symbolic  description  of  the  L  and  U  matrices. 

The  numerical  L  U  factorization  is  performed  by  subroutine  VMNP 
(using  the  code  generated  by  VMSP)  .  Subroutine  VMBP  then  solves  the 


equations  by  forward  and  back  substitution  steps.  Upon  return  from 
VMBP,  the  resultant  solution  vector  is  stored  in  the  original 
right-hand  side  vector  B  of  .the  equations  A  x  =  B.  For  further 
discussion,  reference  [62]  is  recommended. 


CHAPTER  III 


FORMULATION  OF  THE  OPTIMAL  DESIGN  PROBLEM 
AND  SENSITIVITY  ANALYSIS 


3.1  Introduction 

For  constrained  dynamic  systems  with  just  a  few  degrees  of 
freedom,  the  treatment  of  the  optimal  design  problem  given  in 
reference  1,  in  conjunction  with  the  Newtonian  equations  of 

motion  may  be  quite  sufficient.  However,  for  large  systems, 
owing  to  nonlinearity  in  the  system  equations,  the  Newtonian 
approach  is  not  convenient.  On  the  other  hand,  it  will  be  evi¬ 
dent  from  the  next  sections  that  the  Lagrangian  sparse-matrix 
formulation  of  the  equations  of  motion  and  implementation  of  a 
STIFF  (GEAR)  integration  algorithm  are  extremely  advantageous. 

In  the  following,  discussions  will  be  confined  to  two  dimen¬ 
sional  systems.  However , their  extension  to  three-dimensional 
systems  is  quite  straightforward,  at  least  theoretically. 

3.2  Formulation  of  the  Optimal  Design  Problem 

The  optimal  design  problem  for  a  mechanical  system  of  n  bodies 
m  joints,  and  s  spring-dampers  in  the  notation  of  Chapter  2  can  now 
be  stated  as  follows:  Determine  be  RP,  p  being  a  positive  integer, 


to  minimize  the  cost  functional 
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^0  =  gQ(b,z(0) ,1(0),  z(t1) , i(t^)) 


J0 


Lg[t,z(t),y(t),£(t),b]dt 


(3.1) 


where 


&  [ (4 2  9  J  0jl>***j(s“X)] 

T 

[&]_> &5>  •  •  •  >*-(48  _  3)  ] 


(3.2) 


subject  to  the  following  conditions: 

(a)  state  equations  of  motion  of  Eq.  (2.33): 

F(t,z,z,£,p,b)  =  P(b)z  +  f (t,z,£,g,b)  =  0  ; 

(b)  constraints  of  Eq.  (2.35): 

$  (z,b)  =  0  ,  a=l,2, . . . ,2m  ;  (3.3) 

a 

(c)  the  equations  related  to  spring-damper  pairs  (Eq.  (2.36)); 

(d)  initial  conditions 
0(z(O) 


r  z< 

,«-(0),b)  5  { 

U< 


'z(O)  -  v(b)  =  0 

;(o)  -  v(b)  =  o 

(e)  functional  equality  constraints 
=  gg(b,z(0)  ,1(0)  ,z(t1)  ,I(t1)) 

+ 


(3.4) 

(3.5) 


rcl 


'0 


L  [t,z(t)  ,p(t)  ,£,(t)  ,b]dt  =  0  ,  8=1, ...,r'  ; 

P 

(3.6) 


and/or  functional  inequality  constraints 

=  gg(b,z(0),£(0)  .z^)  ,&(t1) ) 

fzl 

L  [t,z(t),jj(t),n(t),b]dt  <_  0  ,  8=r'+l,...,r  ; 

P 

(3.7) 


4 


* 


(f)  design  parameter  constraints 
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s* =  If sb  +  rfw  S2<0)  +  5!m  “(0)  +  5z(ti> 


^Tsi(ti)+[0  [tl  Sz+f7511 


3L.  .  3L  ,  {  , 

"di61  +Fb5b  dt 


(3.13) 


The  first  variations  of  Eqs.  (3.4)  and  (3.5)  give: 


6z(0)  -  5b  =  0 


«A(0)  -  —  fib  -  0 


(3.14) 


To  express  6\p  solely  in  terms  of  6b,  an  adjoint  variable 
X(t)  =  [X^,X^, X*  is  introduced  through  the  identities  [1,70] 
obtained  from  Eqs.  (2.33),  (2.35)  and  (2.36): 


X  [P(b)z  +  f(t,z,y,£,b)]  *  0 


(3.15) 


X^  $>(z,b)  =  0 


(3.16) 


X,T[I  £  +  C(z,£,b)]  =  0 


(3.17) 


where 


X  -  (X1,X2» • • • >x6n) 

a ess  ' 

X  =  (X6n+l’**',X6n+2m) 


X'  (X6n+2m-KL,*"’X6n+2m+4s) 


(3.18) 


Integrating  Eqs.  (3.15),  (3.16),  and  (3.17)  from  0  to  t^,  one 


obtains 


fa  _T 

X  [P(b)z + f (t,z,£,p,b) ]dt  =  0 

Jo 


*3-  X^  $(z,b)dt  =  0 
J0 


(3.19) 


(3.20) 
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C1  X  ,T[I  £  +  C(z,4,b)]dt  =  0 
Jn 


(3.21) 


Integrating  the  first  term  in  Eq.  (3.19)  by  parts,  one  gets 


XTP(b)z 


t.  (t-. 

-  [iT  P(b)z-XT  f(t,z,A,li,b)]dt  =  0 

0  J0 


(3.22) 


Similarly,  integrating  the  first  term  in  Eq.  (3.21)  by  parts,  one 


X'T  II 


|fcl  fl 

[j 

Jn 


>,T  ll-  X,T  Z(z,z,  b)]dt  =  0  (3.23) 


'0  ;0 


The  first  variations  of  Eqs.  (3.22),  (3.20),  and  (3.23)  lead  to 

^ASablal4b+,(b)Sl]|*i  .  J^ipT(i(j^lsb 
+  p(bH  -xTif62-xT^6,-  xTf 


-  *T  -f  Sbldl:  - 0 


x  {z+lt 6b]  dt  ■ 0 

:[I  ft]  1  -  r  i,T(X  6  I)  -  X’ 

o  Jo  L 

-  x’TH6bJdt  ■  ° 


(3.24) 


(3.25) 


T  35  .  ,  ,T  35 

al6z  -  X 

(3.26) 


Integrating  the  first  term  under  the  integral  in  Eq.  (3.24) 
by  parts  and  making  necessary  adjustments,  one  has 


XT  P(b)6z 


t  -  f  F 


*1  P-T„,,  N  .  rT  3  (P(b)z)  ,, 
X  P(b)6z  -  X  3  b  -  6b 


-T  3f  .  rTDf  .  rT  3f  ..  rT 

d76z_x  317  6,1  '  x  3^  "  x 


Hb_ 


dt  =  0  (3.27) 


It  is  to  be  noted  here  that  in  the  second  term  under  the  integral 
of  Eq.  (3.27)  the  differentiation  is  of  P(b)  only  with  respect  to 


b  and  -f  should  not  be  substituted  for  P(b)z,  i.e., 
if  P(b)z  =  P_j_(b)£j  ,  j  summed. 


3 (P(b)z)  „  „  .  „ 

~3b -  5b  =  Pij,kZj  5bk  *  J>  k  summed. 


where 


*Vb) 


ij  ,k  3  b, 


(3.28) 


It  is  further  noted  that  6z(0)  and  6£(0)  are  expressed  in  terms 
of  6b  by  the  Eq.  (3.14). 

Now  to  express  the  4th  and  5th  terms  of  the  right-hand  side  of 
Eq.  (3.13)  in  terms  of  6b,  the  following  boundary  conditions  are 


introduced : 


(3.29) 


X2+4j  =  (tx)  ’  j  0.1.2,. ..,(s-l)  J 

1+4  j 

Then,  with  Eqs.  (3.14)  and  (3.29),  Eqs.  (3.27)  and  (3.26)  can  be 


rewritten  as 


XT(0)  PCb)  f£  6b  -  JbCtp  +{oX[xT  P(b)6 


-  **  H«  - 


§H 


dt  =  0 


(3.30) 


Sr1  X'(0)  61  <0)  =  SyX  X'(0)3-J±i6b  =  SyX  6*  (ti> 

X2+4J  «U*J  -  x2+4j  3b  l0  (t,)  "l-MJ 

J  _  J  J  1+4  .? 


-rt- 


T(I6£)  -  X,T  6z  -  X’T  5 £  -  X’ 

3  2  a  36 


*  If »} 

(3.31) 
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Now  let  X,  X  ,  and  X'  satisfy  the  adjoint  equations 

F*  =  -PT  +  Ml  x  +  ^X  +  X’  -  ^  -  0  (3.32) 

dt  3z  3z  3z  3z 


T  T 

,.s 


3y 


and 


T  T  T 

_  -T  dX'  ,  3f  T  ,  ar  v  9L_ 

^  -  1  'dT  +  _3r  X  +  3£  X  3* 


=  0 


(3.33) 


(3.34) 


with  initial  conditions  (3.29)  already  chosen.  With  Eqs.  (3.32), 
(3.33),  and  (3.34),  one  obtains  from  Eq.  (3.30) 

xT(0)  P(b)  f  +{o1[frff^  +  x’TflS2 


3  L 
3z 


Sz-~Su-  i'T  m  +  X'T-?|s£-S 


3y 


3  £ 


3  £ 


_  3(P(b)z)_  6b  _  XT  |i  fibldt  =  0 

8  b  o  b 


(3.35) 


Hence  from  Eqs.  (3.25),  (3.31),  and  (3.35),  one  obtains, 

s— 1  3v. 


xT 


(0)  P(b)  6b  +  X.W)  3z(f1)"  6z(tl) 


s-1 


-  I 


JLs_ 


5£ 


(ti) 


**»&!  1+4J 


rc(- 

0  L\ 


3  L  .  3  L  j,  3  L 

57 5z  Sl 


-T^Sb  Sb  6b 


3b 


-X"T 


§H 


3b 
dt  =  0 


3b 


(3.36) 


Therefore, 


3z(t1) 


6z(tl)  +T^^6£(t1)  + 


3T(t1) 
3  v 


3  L  .  .  3  L  ~  ,  3  L  ~ .  m 

—  Jz+—  Sp+jjSl  |dt 


S_1  ,  (0)  8v.i+l 

3b 


=  r(0)  P(b)  f£«b  +  tIQ  ^  ' 2+4 j 

f1! 


6b 


X-T  6b  +  X-T  II  ab+r1  II  6b  +X  ,T  If  6b  kt 


o 


3b 


3b 


3b 


3  b 
(3.37) 
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Eq.  (3.37)  can  now  be  used  to  eliminate  explicit  dependence 


of  Sip  on  6  z,  6y,  and  6£  and  the  steepest  descent  programming  meth¬ 
od  can  be  used  to  carry  out  sensitivity  analysis  and  iterative 
optimization.  Substituting  from  Eq.  (3.37)  into  Eq.  (3.13),  one 


obtains 


jii k  +  PCb). 


‘lb5b+Lrlw  +  'lT(0)  PCb)]f  5b 
+  stY-i^-+x'(o)p±i8b+ftir^- 


jT  3(p(b)z) 


(3.38) 


As  stated  previously,  Eq.  (3.38)  can  be  applied  to  all  ip  , 

p 


8=0,1, ... ,r.  For  convenience,  define 


T  3  %  +  T  3  %  V(0)  Pfb)  1  is  +  °fT  38e 
,  +  L3z<°)  <W  J  3b  AJ.3J,  (0) 

—  J  l-f-A  j 


iV  ,13u 


s-lr  3g, 


+  >'  p 

+  X2+4j 


(3.39) 


which  is  the  design  sensitivity  coefficient  vector  of  iftg  with  res¬ 
pect  to  the  design  parameter  b.  With  Eq.  (3.39),  one  can  write 


=  o*0 


6b 


(3.40) 


Sip0  =  6b  ,  8=1,2, ...  ,r  . 

P 


(3.41) 


Eqs.  (3.40)  and  (3.41)  now  provide  6i^n  and  $\p  solely  in  terms  of 

U  8 

fib.  The  generalized  steepest  descent  or  gradient  projection  method 


of  [1]  may  now  be  applied  to  carry  out  constrained  sensitivity 
analysis  and  iterative  optimal  design. 


3.4  Comparison  of  the  Corrector  Equations  for  the 
Equations  of  Motion  and  the  Adjoint  Equations 

For  system  Eqs.  (2.33),  (2,35),  and  (2.36),  the  corrector 

equation  can  be  written  (see  Chapter  II)  as: 


_L_  p(ii)  +11  11 

a  h  K  J  3z  3y 
o 


M.  0 

3  z 


3z 


0 


3l 

3 1 

Az 

~  — 

F 

0 

AU 

=  - 

$ 

-i-i+M. 
eQh  AT3Jl 

LZ 

I 

(3.42) 


In  a  similar  manner,  the 
Eqs.  (3.32),  (3.33),  and 


_i_pT  +  3f_ 

80h  +  3z 

M1 

3u 


corrector  equations  for  the  adjoint 
(3.34)  can  be  written  as: 


r  \T 

/  \T 

“  “ 

I  —  ) 

(1L) 

AA 

F* 

[dzj 

\32/ 

0 

o 

AA 

*  - 

0 

_L_iT+llT 

SQh1  3JL 

aa' 

V 

(3.43) 


where  the  following  relations  have  been  used: 
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9  F  _  3  f 
3  z  3  z 


9_F  _  9_f 
34  9  4 


9  F  _  9  f 
9p  9  vi 


(3.44) 

'From  Eqs.  (3.42)  and  (3.43)  it  is  clear  that  the  coefficient 


matrix  of  Eq.  (3.43)  is  the  transpose  of  the  coefficient  matrix  of 


Eq.  (3.42),  except  for  the  minus  sign  before 


SQh 


in  the  diagonal 


terms  of  Eq.  (3.42).  However,  since  the  adjoint  equations  are  inte¬ 
grated  backwards  in  time,  the  negative  time-step  will  make  the 
second  matrix  exactly  the  transpose  of  the  first.  Thus,  essential 
computational  advantages  accrue. 


3.5  The  Solution  of  the  Adjoint  Equations 


The  system  of  adjoint  equations  (3.32),  (3.33),  and  (3.34) 
are  to  be  solved  with  the  terminal  conditions  (3.29).  With  the 
substitution 


t'  =  t,  -  t 


(3.45) 


the  adjoint  equations  become 


V  =  pTdf,  *-  +  hI-  +  mI 

1  dt  9z  3z  9  z 


A'  - 


9  L 
9  z 


T  T 

4  .  s  Lf_x-2£--o 

1  9y  9y 

l  •  =  fT  dT  3f^  -  3^  ,  _  3l!  =  o 
41  "  dt1^  34  A  94  X  94 


0  (3.46) 

(3.47) 

(3.48) 
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\ 


with  the  initial  conditions: 


rI<0)  -  rfo-)  p(brl 


X’(0) 
2+4  j 


j=o,i. 


(8-1)  | 


J 


(3.49) 


Then  the  corrector  coefficient  matrix  becomes  exactly  the  trans¬ 
pose  of  the  original  coefficient  matrix.  To  make  full  use  of  the 
results  obtained  in  the  process  of  solution  of  the  original  set 
of  equations  of  motion,  the  time  instants  (indexed) ,  Gear  constant 
3^,  step-size,  the  solution  variables,  and  the  coefficient  matrix 
elements  are  stored  in  a  direct  access  disk  at  each  time  grid. 

The  matrix  elements  or  the  G-vector  of  ADAMS  2-D  [14]  are  origin¬ 
ally  arranged  with  column-wise  representation. 

To  make  use  of  the  same  subroutine  as  in  the  SPARSE-MATRIX 
package  for  LU  factorization  of  the  transposed  coefficient  matrix, 
some  modifications  have  been  made  in  the  subroutine  S08000  of  ADAMS 
2-D  [14]  to  represent  the  transposed  matrix  in  a  column-wise  manner. 
It  is  then  stored  in  the  direct  access  device.  It  is  to  be  noted, 
however,  that  the  time  grid  of  the  backward  integration  for  the  sol¬ 
ution  of  the  adjoint  equations  will  not,  in  general,  coincide  with 
the  original  grid.  So  interpolation  of  the  original  solution  var¬ 
iables  and  the  second  components  of  the  Nordsieck  vector  [6,2]  is  re¬ 
quired  for  calculation  of  the  right-hand  side  of  the  corrector  equa¬ 
tions  of  the  adjoint  set  and  the  integrands  of  the  sensitivity 
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matrices.  This  is  done  by  a  subroutine  INTERP  built  in  the  sub¬ 
routine  DYNANL  of  ADAMS  2-D  (also  see  Chapter  V) .  SOINEW  is  the 
subroutine  for  calculation  of  the  right-hand  side  of  the  new  cor¬ 
rector  equations,  which  is  built  into  S01000  of  ADAMS  2-D. 

The  G-vector,  however,  is  not  interpolated.  It  is  approx¬ 
imated  by  its  value  at  the  nearest  original  time  grid  point.  It 
can  actually  be  kept  unaltered  for  several  small  time  steps;  in 
particular,  during  the  CNTRLT  and  Newton  iteration  operations  in 
DIFSUB  of  ADAMS  2-D.  Finally,  ADJONT  is  a  version  of  ADAMS  2-D 
that  solves  the  adjoint  equations,  with  all  above  considerations 
built  in.  Chapter  V  deals  with  the  subroutines  in  more  detail. 

3.6  Static  Sensitivity  Analysis 
for  the  Solution  Variables 

From  Eqs.  (3.39),  (3.40),  and  (3.41)  it  is  observed  that  for 
the  calculation  of  sensitivity  coefficients  one  needs  derivatives 
of  the  initial  values  of  some  solution  variables  with  respect  to 
the  design  parameters.  These  are  obtained  as  follows. 

For  static  equilibrium  the  system  Eqs.  (2.33),  (2.35),  and 
(2.36)  reduce  to  the  following  set  of  algebraic  equations: 

f(z,i,,y,b)  =  0  (3.50) 

$(z,b)  =  0  (3.51) 

C(z,il,b)  =  0  (3.52) 

Differentiating  these  equations  with  respect  to  b  one  has 

9z  9b  +  9p  9b  +  9 1  9b  ~  9b 


(3.53) 
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3£  3_z 
3z  3b 


3$ 

ab 


H  3z  ac  jn  _  h 

dz  3b  35-  3b  “  3b 


(3.54) 

(3.55) 


Eqs.  (3.53),  (3.54),  and  (3.55)  represent  a  set  of  (3n+2m+4s) 

|p  |p  and  ||  .  In  matrix  form  these  eq- 


linear  equations  for 
uations  are  equivalent  to 


3v  3F 
J1  3b  =  3b 


(3.56) 


where  (using  Eq 

.  (3.42)), 

"3f 

3f 

3f“ 

3z 

3y 

3 1 

Jis 

_3$ 

3z 

0 

0 

(3.57) 

H 

0 

H 

_3z 

3  Z_ 

is  the  Jacobian  matrix  of 

the  static  analysis  and 

3v  _ 
3b  = 

r3z 
_3b  ’ 

3b 

3  Jl"|T 

’  3bJ 

(3.58) 

3F  = 
3b  " 

r  3f 

3  b 

i  “ 

3$ 

3b  * 

~  3  bj 

(3.59) 

The  sparse  matrix  codes  generated  during  static  analysis  are 
utilized  to  solve  for  ||  ,  with  the  help  of  the  subroutine  DIFSUB. 
The  static  sensitivity  analysis  results  for  a  spring-reset  plow¬ 
share  mechanism  (see  Chapter  VI)  are  given  in  Tables  6.7  and  6.8 


in  Chapter  VI. 


CHAPTER  IV 


OPTIMAL  DESIGN  ALGORITHM 


4.1  Steepest  Descent  Method  with  Constraint 
Error  Compensation 

The  general  optimal  design  problem  formulated  in  Chapter  III 
can  be  solved  by  the  generalized  steepest  descent  programming 
technique  presented  in  Ref.  [1].  Here  the  technique  will  only 
be  discussed  very  briefly. 

After  the  sensitivity  analysis  of  Chapter  III,  the  problem 
is  reduced  to  finding  a  vector  6b  that  minimizes  and  corrects 
constraint  violations.  For  this  purpose,  the  following  defini¬ 
tions  are  made.  Define  a  set  of  indices 


A  =  {S|ij>s  +  e  >_  0} 


(4.1) 


and  a  column  vector  of  e-active  elements  of  ipg  , 


4>  = 


Se  A 


(4.2) 


In  order  to  assure  the  satisfaction  of  the  constraints,  the  viol¬ 
ations  will  be  corrected  by  demanding  that 


6ij>  _<  Aty  g 


(4.3) 


where  SeA  and  Aijig  is  a  desired  change  in  the  constraint  function 
\J»g.  Generally,  is  chosen  to  be  -i|>g. 

In  the  notations  of  Chapter  III,  the  problem  reduces  to 


finding  5b  to  minimize 

*oT 

=  Ji  u  5b 

subject  to  the  constraints 

TT 

5$  =  J?/  5  b  _<  Arp 

and  the  quadratic  step  size  constraint 
5bT  W  5b  _<  n2 


(4.4) 


(4.5) 


(4.6) 


where  n  is  a  small  number  and  W  is  a  positive  definite  weighting 
matrix. 

The  Kuhn-Tucker  necessary  conditions  of  nonlinear  programming 

may  now  be  applied  to  solve  this  reduced  problem.  According  to 

the  Kuhn-Tucker  Theorem,  it  is  necessary  that  there  exist  a  scalar 

multiplier  y  >  0  and  a  vector  of  multipliers,  y  =  { yc ,  S e  A}  , 

U  b 

Yg  >_  0  for  inequality  constraints,  such  that 
~T 

yT(Jl^  5b  -  A$)  =  0  (4.7) 


and 


3  G 

3  (5b) 


=  0 


(4.8) 


where 


T  ~T 

G  =  5b  +  yT(^  <Sb  -  Aiji)  +  yQ(5bTW5b  -  n2)  (4.9) 


Equations  (4.8)  and  (4.9)  give 


Solving  for  6b,  one  gets  from  Eq.  (4.10) 


corresponding  constraint  6ij>g  -  Aip s  _<  0  should  have  been  strictly 
satisfied.  Therefore,  should  be  removed  from  4! .  Then  a  new 
reduced  set  ij;  is  formed  and  the  multipliers  of  this  new  set  are  re¬ 
calculated  . 

From  Eqs.  (4.17)  and  (4.11)  one  obtains 

fib  -  -  +  y1)  -  W-1  $  y2  (4.18) 

Y0 

Defining 

fib1  =  +  y1)  (4.19) 

fib2  =  -  W-1  $  y2  (4.20) 

the  change  fib  of  Eq.  (4.18)  can  be  written  as 

6b  *  -  6b1  +  6b2  (4.21) 

^0 

In  practice,  instead  of  choosing  q  in  Eq.  (4.6),  Yq  is  chosen 
directly  to  give  the  step-size  . 

The  following  relations  can  be  easily  shown  to  follow  as  a 
necessary  consequence  from  the  above  equations: 


±  2 
fib  Wfib  =  0 

(4.22) 

2 

V  fib  =  Aip 

(4.23) 

~T 

$  fib1  =  0 

(4.24) 

-  fib1  _<  0 

(4.25) 

From  these  relations  it  can  be  observed  that: 
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section  gives  the  following  optimal  design  algorithm: 

(j) 

Step  1:  Estimate  the  optimum  design  parameter  vector  b 

and  solve  the  state  equations  of  motion  (2.33),  (2.35),  and 
(2.36)  for  z^(t),  V^(t),and  A^\t)  corresponding  to  b  ^  with 
SPARSE  MATRIX  and  STIFF  INTEGRATION  (GEAR)  algorithms  imple¬ 
mented  through  the  ADAMS  Program  and  at  each  time  step  store 
the  solution  variables,  indexed  time-instants,  time  step-size, 
and  corrector  coefficient  matrix  in  a  direct  access  disk. 


Step  2;  Check  constraints  of  Eqs.  (3.7),  (3.8),  and  (3.10) 
and  form  the  vector  of  constraint  functions  if  ,  consisting  of 
such  that  >_  -  e,  where  £  >  0  is  small.  Also  choose  the 
the  constraint  error  correction  vector  A  if  =  -  if  . 

Step  3:  Implement  a  modified  ADAMS  program  to  solve  the  ad¬ 

joint  differential  equations  (3.46),  (3.47),  and  (3.48)  with 

$  rv  Ip  o 

initial  conditions  (3.49),  to  obtain  A  u(t)  and  A  p(t),  ^ 

p 

corresponding  to  e -active  constraint  functions. 


Step  4 :  Compute 


,*0  3g0  [~3§0  r^0Tf0srasl9v  Sy1  1"  3 g0 

1  -tt  +  l^(o)-+x  ((Wb)Jrb  +  ji0L^7W 


+  A'2+Aj(o) 


Jo  L3b  - 


I^0T  3(P(b)z)  _  r^0T  if 


=^0T  ii  ,,  ^  35 
"  A  9b  "  x  9b  ‘ 


(4.27) 


3b 


3i!L  +  xV(omb)]f^+T 


»S, 


j=0  '-dz 


(0) 

l+4j 


+  X 


*rT 


*  ^3 


3v . 

An 

-  J 
3b 

+ 

Jo  - 

f  !I6  TV  3(P(b)  &) 

JL  Lab  A  3b 


if  f^gT  ii  x;  il 

A  3b  "  A  3b  ~  A  3b_ 


dt 


(4.28) 


and 


M 


# 


|Jl^  W  ^  Z^ ,  if  ip  is  not  empty 


,  if  »P  is  empty 


(4.29) 


jh  ’Pa 

where  Z  is  the  (Pxg)  matrix  constituted  of  the  vectors  Z  p, 

p  is  the  number  of  design  parameters  and  8  is  the  number  of 

violated  constraints,  and  W  is  a  weighting  matrix  [1] . 

Step  5:  Choose  step-size  “  >  0  and  calculate  the  Lagrange 


0 


multiplier  vector  y  from 

M  Y1  =  -  a*  w*"1  /° 
’p’p 


m-A* 


(4.30) 

(4.31) 


and 


1  .  -  2 
Y  =  Y  +  2yqY 


(4.32) 


Step  6:  Check  the  algebraic  sign  of  each  component  of  y 

associated  with  inequality  constraints.  If  any  components 
of  y  are  negative,  redefine  ip  by  removing  the  corresponding 
terms  from  \p  and  return  to  Step  4. 


75 


Step 

1  2 

Compute  6b  ,  6b  ,  and  6b  from 

6b1 

=  w"1(/°+^y1) 

(4.33) 

6b2 

=  -  w_1  Jt'V 

(4.34) 

and 

5b  = 

»  -  6b1  +  6b2 

2Yq 

(4.35) 

5b  gives  a  design  improvement. 

Step  8 :  Compute  b^+1^  =  b^  +  6b  (4.36) 

Step  9:  If  the  constraints  are  satisfied  and  ||6b^||  is 

sufficiently  small,  terminate  the  process.  Otherwise,  re¬ 
turn  to  Step  1  with  b^+1^  as  the  best  estimate  of  the  design. 
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CHAPTER  V 

THE  DYNAMIC  ANALYSIS  AND  DESIGN  SYSTEM  (DADS)  PROGRAM 


5.1  Introduction 

The  Dynamic  Analysis  and  Design  System  (DADS)  computer  program 
is  developed  to  carry  out  the  dynamic  analysis,  design  sensitivity 
analysis,  and  optimal  design  formulations  described  in  Chapters  II, 
III,  and  IV  for  general  nonlinear  mechanical  systems.  Provisions 
for  regenerating  sparse-matrix  codes  at  necessary  time  instants  of 
dynamic  and  adjoint  analyses  have  been  made  so  that  the  program  can 
handle  systems  with  intermittent  motions  (see  Plow-share  mechanism 
in  Chapter  VI  (Section  6.3))  with  sufficient  ease. 

The  DADS  program  executes  in  two  main  phases:  (1)  The  dynamic 
analysis  phase  and  (2)  The  sensitivity  analysis  and  optimization 
phase.  The  dynamic  analysis  phase  of  the  program  generates  the 
sparse  matrix  code  for  pivoting  and  Lu  factorization  of  the  Jacobian 
matrix  and  solves  the  differential-algebraic  equations  for  the  state 
variables  during  a  specified  time  interval.  It  employs  sparse 
matrix  codes  and  the  numerical  integration  routine  presented  in 
Chapters  II  and  III.  It  then  stores  the  Jacobian  matrix  and  state 
variables  on  a  direct  access  disk  for  use  in  the  sensitivity  analysis 
phase.  The  dynamic  analysis  phase  also  performs  static  sensitivity 
analysis  for  the  solution  variables  (see  Section  3.6). 

The  sensitivity  analysis  and  optimization  phase  of  the  program 
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incorporates  the  techniques  of  Chapter  IV  into  various  subprograms 
and  solves  adjoint  equations  (see  Section  3.3)  using  the  data  stored 
on  disk.  The  same  numerical  integration  subroutine  is  used  to 
solve  the  adjoint  equations.  The  program  further  computes  sensiti¬ 
vity  coefficients  and  the  necessary  design  improvements.  The  pro¬ 
cess  is  repeated  until  optimum  results  are  achieved. 

5.2  Main  Features  of  DADS  Computer  Program 
The  DADS  program  has  the  following  features: 

(1)  The  mechanical  systems  treated  are  discrete,  nonlinear,  and 
constrainted  (two-dimensional  at  present) ; 

(2)  General  nodal  formulation  of  equations  of  motion  for  the 
bodies  [7]; 

(3)  Necessary  data  being  given,  formulation  of  the  equations  of 
motion  and  the  Jacobian  matrix  is  automatic; 

(4)  Integration  of  a  combined  set  of  differential  and  algebraic 
equations  (DAE)  is  performed; 

(5)  The  following  three  algorithms  are  used: 

(a)  Newton  iteration; 

(b)  Sparse  Matrix  techniques  for  the  solution  of  the  linearized 
simultaneous  equations; 

(c)  Stiff  Integration  Algorithm  (GEAR)  for  numerical  integration. 

(6)  Types  of  analyses  performed : 

(a)  Static  analysis; 

(b)  Transient  analysis; 
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(c)  Static  sensitivity  analysis  for  solution  variables; 

(d)  General  design  sensitivity  analysis; 

(e)  Optimal  design. 

The  static  and  transient  analysis  part  of  the  program  was 
mainly  developed  by  Orlandea  [7]  and  Wehage  [14].  The  correspond¬ 
ing  references  are  recommended  for  further  details.  Fig.  5.1  shows 
an  outline  of  DADS  program  capabilities,  Fig.  5.2  shows  the  main 
subroutines  used  in  the  dynamic  analysis  phase  of  DADS,  and  Fig.  5.3 
identifies  the  main  subprograms  used  in  the  design  sensitivity  anal¬ 
ysis  and  optimization  phase,  and  Fig.  5.4  is  a  flow  diagram  of  the 
DADS  computer  program. 

5.3  Brief  Description  of  the  Dynamic  Analysis  Phase 
The  dynamic  analysis  phase  (DYNANL)  of  the  DADS  program  in¬ 
volves  establishing  the  sparse  matrix  code  description  of  the  mechan¬ 
ical  system  and  solving  the  differential  and  algebraic  equations  for 
the  state  variables.  As  shown  in  Fig.  5.2,  this  involves  two  major 
steps:  (i)  generation  of  an  initial  sparse  matrix  code  (including 
pivoting  and  LU  factorization  code)  and  (ii)  repetitive  solution  of 
linearized  equations  for  the  state  variables  during  the  time  interval 
of  interest. 

In  the  first  step  of  dynamic  analysis,  estimates  of  the  initial 
configuration  of  the  system  are  provided  by  INDATA  and  are  used  by 
VARSET  to  initialize  a  state  variable  vector  for  subsequent  use  by 
the  numerical  integration  routine  DIFSUB.  A  compact  numbering 


DYNANL 


Dynamic  Analysis  Phase 
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Figure  5.4. 


Flow  Diagram  of  the  DADS  Computer 
Program. 
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system  identifying  bodies,  joints,  and  spring-damper  elements  is 

used  to  input  data  through  VARSET  and  provides  the  necessary  des- 

\ 

cription  of  the  mechanical  system  configuration.  These  data  are 
used  to  construct  the  Newton  corrector  state  equation  (3.42)  and  the 
system  of  adjoint  equations  (3.43)  (for  use  in  the  sensitivity  anal¬ 
ysis  phase) .  This  information  is  used  in  S03000  to  generate  init¬ 
ial  vectors  of  row  and  column  indices  for  nonzero  entries  in  the 
Jacobian  matrix  and  to  assemble  the  standard  equations  that  are  re¬ 
quired  on  the  right  hand  side  of  Eq.  (3.42).  Similarly,  descrip¬ 
tions  of  user-supplied  row-column  positions  of  additional  non¬ 
standard  elements  are  provided  by  incorporating  the  necessary  code 
in  USET.  A  symbolic  description  of  the  resultant  matrix  is 
printed  by  DEBUGG  for  reference  purposes  and  a  column-ordering  perm¬ 
utation  vector  is  generated  in  S08000. 

Subroutine  S01000  evaluates  the  Jacobian  matrix  and  right  hand 
side  of  Eq.  (3.42).  Its  purpose  is  to  (i)  evaluate  force  and  dis¬ 
placement  functions  of  time  that  are  provided  by  the  user  through 
subroutine  FOREXT,  (ii)  transfer  the  state  variables  from  a  single 
vector  used  by  the  numerical  integration  routine  to  the  standard  var¬ 
iables  (and  user-supplied  variables  through  U S0LV1) ,  (iii)  evaluate 
that  part  of  the  Jacobian  matrix  that  is  associated  with  the  standard 
equations  (and  user-supplied  equations,  US01V2 ) ,  using  updated  var¬ 
iables  from  step  (ii)  and  the  previously  generated  column-ordering 
permutation  vector  (S08000) ,  and  (iv)  evaluate  the  standard  equations 
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(and  user-supplied  equations,  US0LV3).  Subroutines  LSPTRS,  LSSLV1, 

LSSLV2,  and  LSSLV3  are  included  as  user-supplied  routines  to  incor-  ^ 

porate  nonstandard  highly  nonlinear  springs  and  dampers,  as  described 

in  the  plow  share  mechanism  example  of  Chapter  VI.  Finally  a  sparse 

LU  factored  description  of  the  matrix  is  generated  in  subroutine 

INVERT. 

The  second  step  of  dynamic  analysis  is  then  to  numerically  inte¬ 
grate  the  system  of  equations  during  the  time  interval  of  interest. 

This  is  accomplished  by  the  numerical  integration  subroutine  DIFSUB, 
which  repeatedly  calls  S01000  to  update  the  Jacobian  matrix  and  system 
of  equations  as  it  executes  the  iterative  corrector  formula  of  Eq. 

(3.42).  The  Jacobian  matrix  and  the  results  of  numerical  integration 
are  stored  in  a  direct  access  disk. 

In  sensitivity  analysis,  DYNANL  is  called  in  ADJONT  (adjoint  anal¬ 
ysis  phase)  (Fig.  5.4),  which  is  again  called  in  AHLPSI,  where  the 
sensitivity  coefficients  are  calculated.  For  this  case,  DYNANL  reads 
stored  data  from  the  disk  and  passes  through  two  distinct  steps.  In 
the  first  step,  the  sparse  matrix  codes  for  pivoting  and  LU  factoriz¬ 
ation  are  generated  for  the  transpose  of  the  Jacobian  matrix,  through 
INVERT.  In  the  second  step,  the  system  of  adjoint  equations  (3.32) 
to  (3.34)  are  numerically  integrated  through  DIFSUB,  which  repeatedly 
reads  data  from  the  disk,  calls  S01NEW  for  the  evaluation  of  the 
right-hand  sides  of  the  adjoint  system,  solves  the  adjoint  variables 
at  various  time  steps,  and  stores  the  solution  results  on  a  disk  for 
use  in  the  calculation  of  sensitivity  coefficients. 
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5.4  Description  of  the  DADS  Program 

As  noted  in  Section  5.1,  the  computer  program  DADS  is  consti¬ 
tuted  of  two  main  phases:  (i)  Dynamic  Analysis  phase,  and,  (ii) 
design  sensitivity  analysis  and  optimization  phase.  The  main  flow 
of  the  program  is  described  below  (see  also  Fig.  5.4). 

In  the  main  program,  the  initial  estimates  and  bounds  on  de¬ 
sign  parameters,  number  of  constraints,  stepsize,  percentage  of 
cost  function  reduction,  and  other  parameters  related  to  the  design 
problem  (see  Chapters  III  and  IV)  are  read.  Then  system  parameters, 
such  as  masses,  moments  of  inertia,  locations  of  centers  of  masses, 
applied  constant  forces,  joint  types,  and  spring-damper  parameters 
are  read,  through  the  subroutine  INDATA.  The  subroutine  RELATE 
relates  the  variables  and  parameters  of  the  dynamic  analysis  (ADAMS  2 
and  DYNANL;  see  also  Section  2.1  of  Chapter  II)  to  the  updated  de¬ 
sign  parameters  (Chapters  III  and  IV). 

For  dynamic  analysis ,  DYNANL  is  called  through  ADAMS2.  The 
Jacobian  matrix,  solution  variables,  time  step,  time  instant,  and 
order  of  numerical  integration  are  then  stored  in  a  direct  adcess 
disk.  They  are  first  used  to  evaluate  the  integrands  of  the  cost 
and  constraint  functionals  of  Eqs.  (3.1),  (3.6),  (3.7),  and  (3.11), 
through  subroutine  HALS. 

Next,  the  inequality  functional  constraints  of  Eqs.  (3.7)  and 
(3.10)  are  tested  and  the  corresponding  adjoint  equations  (3.32)  to 
(3.34)  are  solved  through  INFUNC,  TEQFUN,  AHLPSI,  and  ADJONT.  At 
this  stage,  DYNANL  (Adjoint  Analysis  Phase)  is  called  by  subroutine 


ADJONT,  The  matrix  of  sensitivity  coefficients  l  (see  Eq.  (3.41)) 

is  calculated  for  the  violated  functional  constraints  through  AHLPSI. 

The  subroutine  DPARMC  then  tests  the  design  parameter  constraints 

and  calculates  the  corresponding  design  sensitivity  vectors.  The  sensi- 

< pn 

tivity  vector  l  u  of  Eq.  (3.40)  is  then  calculated  through  AHLPSI.  The 

subroutine  NEWB  then  computes  the  matrix  M  of  Eq.  (4.29),  solves  Eqs. 

1  2 

(4.30)  and  (4.31)  for  y  and  y  ,  and  computes  the  design  change  6b 
given  by  Eqs.  (4.33)  through  (4.35). 

At  this  stage,  convergence  criteria  are  tested.  If  they  are  sat¬ 
isfied,  the  new  design  is  taken  as  the  optimum  one.  Otherwise,  the 
process  is  repeated  with  the  new  design  paramters  used  as  the  initial  de¬ 
sign  estimate. 

In  the  following  sections,  descriptions  of  the  principal  program 


variables  and  subroutines  are  presented. 


5.4.1  Principal  Variables 
a)  Dynamic  Analysis 

TIMPU  -  current  time 

NB  -  number  of  bodies  including  fixed  body  1 

NJ  -  number  of  joints 

NSD  -  number  of  spring-damper  pairs 

NSDV  -  twice  the  number  of  nonlinear  spring-damper  pairs 

NSDT  -  number  of  torsional  spring-damper  pairs 

EPS  -  maximum  one-step  error  in  numerical  integration 
M  -  vector  of  masses  of  the  bodies 
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JIN  -  vector  of  moments  of  inertia 

X,  Y  -  coordinates  of  the  CM's  with  respect  to  inertial 

reference  frame 

PHI  -  angular  displacements  of  X-axis  of  bodies  with 
respect  to  inertial  refernece  frame  X-axis 

FX,  FY  -  components  of  the  force  applied  at  CM  parallel  to 
the  inertial  reference  frame  axes 

TQ  -  applied  torques 

JT  -  type  of  the  joints 

IB (1,1) ,IB(2,I)  -  numbers  given  to  the  two  neighboring  bodies 
connected  by  the  ifc^  joint 

XI,  Y1  -  co-ordinates  of  the  revolute  joint  or  a  point  on 

the  axis  of  the  translational  joint  with  respect 
to  the  axes  of  the  first  body 

X2,  Y2  -  same  as  above  with  respect  to  the  2n<^  body 

IBSD(1,I) ,IBSD(2,I)  -  numbers  given  to  the  two  bodies  connected 
by  the  ic^  spring-damper  pair 

XF1,YF1  -  co-ordinates  of  the  point  of  attachment  on  the  first 
body  with  respect  to  the  body  fixed  axes 

XF2,YF2  -  same  as  above  with  respect  to  the  2nc^  body 

SK  -  spring  constants 

DC  -  damping  coefficients 

SDL  -  deformed  spring  lengths 

SDLO  -  undeformed  spring  lengths 

SDF  -  constant  forces  applied  along  spring-damper  pairs 

IBSDT(1,I) ,IBSDT(2,I)  -  numbers  given  to  the  two  bodies  con¬ 
nected  by  the  i*-*1  torsional  spring-damper  pairs 

SKT  -  torsional  spring  constants 

DCT  -  torsional  damping  coefficients 
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PHIO  -  initial  difference  of  the  angular  displacements 
of  the  bodies  for  the  undeformed  torsional 
spring-dampers 

TQO  -  constant  couples  applied  at  the  torsional  spring- 
damper  pairs 

V  -  the  time-derivatives  of  SDL’s 

FFX,FFY  -  components  of  the  spring -damper  forces  with  res¬ 
pect  to  inertial  reference  frame 

'  TQS  -  torsional  spring-damper  couples 

UX,UY  -  components  of  the  velocities  of  the  CM's  with  res¬ 
pect  to  inertial  reference  frame 

UP  -  angular  velocities  of  the  bodies  with  respect  to 

inertial  reference  frame 

LM  -  Lagrange  multipliers  of  the  dynamic  analysis 

G  -  vector  of  column-ordered  nonzero  elements  in  the 

Jacobian  matrix 

CL  -  right-hand  sides  of  the  corrector  equations 

JSIZ  -  size  of  the  Jacobian  matrix 

NPOS  -  pointer  to  consecutive  nonzero  positions  in  the 
Jacobian  matrix. 

b)  Sensitivity  Analysis  and  Optimization  phase 

The  principal  variables  occurring  in  the  program,  over  and  above 
those  mentioned  above  are  described  below. 

NB1  -  number  of  design  parameters  (total) 

NAB  -  number  of  artificial  design  parameters 

MS  -  expected  size  of  constraint  set 

NFCE  -  number  of  pure  equality  functional  constraints 

NFCET  -  number  of  transformed  equality  functional  con¬ 

straints 


* 
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NFCI  -  number  of  inequality  functional  constraints 

NDCI  -  number  of  design  parameter  constaints 

NP1  -  number  of  time  grids  in  stiff  integration  in 

ADAMS  2 

B  -  design  parameters 

BI  -  starting  values  of  the  design  parameters 

BL  -  lower  bounds  of  the  design  parameters 

BU  -  upper  bounds  of  the  design  parameters 

COSTJ  -  cost  function 

EPE  -  tolerance  for  equality  functional  constraints 

EPI  -  tolerance  for  inequality  functional  constraints 
EPDI  -  tolerance  for  design  parameter  constraints 

ER  -  convergence  criterion 

IL1M  -  limit  of  iterations  in  optimization  program  step 

W  -  weighted  matrix  (generally  diagonal) 

NP2  -  number  of  time  grids  in  stiff  integration  in 
ADJONT 


AIM  -  adjoint  variable 

GL  -  vector  of  column-ordered  nonzero  elements  in  the 

transposed  Jacobian  matrix 


MPOS  -  a  vector  playing  the  same  role  in  adjoint  anal¬ 
ysis  as  NPOS  in  dynamic  analysis  (see  Chapter  II) 


NPOSRR 

JP1A 


-  a  copy  of  original  NPOSR 

-  a  copy  of  JP1 


(see  Chapter  II) 


ALS  -  the  integrands  of  the  cost  and  constraint  func¬ 
tionals 


GS  -  non-integral  parts  of  the  cost  and  constraint 

functional 
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AHSMT  -  maxima  of  the  integrands  of  the  transformed 
equality  constraints 

DLDZ  -  derivatives  of  the  integrands  of  the  cost  and 
constraints  functionals  with  respect  to  the 
state  variables  (both  primary  and  secondary) 

HT-T  -  sensitivity  coefficients  for  the  cost  functional 

HLFS1  -  sensitivity  coefficients  for  the  constraints 

functionals 

PSIC  -  constraint  functionals 

NCV  -  number  of  constraint  violations 

NCVID  -  indices  of  the  violated  constraints 

Some  other  variables  are  defined  during  the  description  of  the 
subprograms  in  the  next  section. 

5.4.2  Description  of  the  Subprograms 

In  this  subsection  a  brief  description  of  each  of  the  main  sub¬ 
programs  of  DADS  is  given  and  some  of  the  call-list  variables  not 
appearing  in  the  previous  subsection  are  explained  briefly. 

SUBROUTINE  INDATA  (IECH0,  ALM) 

The  subroutine  INDATA  reads  input  data  that  characterizes  the 
mechanical  system.  It  also  reads  initial  values  of  adjoint  variables 
and  the  maximum  values  of  nonlinear  spring-damper  pairs  into  the  pro¬ 
gram. 

IECH0  -  a  logic  variable  for  writing  the  data  on  paper 

■  0,  data  is  not  written, 

■  1,  data  is  written. 

SUBROUTINE  RELATE  (YY,  IMODE) 

The  subroutine  RELATE  relates  the  parameters  and  variables  of  the 


if 
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dynamic  analysis  to  the  updated  design  parameters. 

YY  -  an  array  which  contains  the  dependent  variables  (dynamic 

» 

analysis)  and  their  derivatives.  Its  dimension  should  be  at  least 
(JSIZ,7)  where  JSIZ  is  the  Jacobian  size. 

IMODE  -  a  logic  variable  characterizing  the  type  of  analysis  per¬ 
formed, 

=  0,  neither  static  nor  dynamic  analysis, 

=  1,  static  analysis, 

=  2,  dynamic  analysis. 

SUBROUTINE  ADAMS 2 ( IFLAG1 ,  SOLNO,  JSIZ,  DLDB,  DFDB,  DPHDB , 
DZIDB,  YY,  PB,  ALM,  DNDB) 


The  subroutine  ADAMS2  feeds  in  some  flags,  time  step,  maximum 
and  minimum  time  steps,  and  calls  subroutine  DYNANL  to  initiate  tran- 
ient  solution. 

IFLAG1  -  a  logic  variable, 

=  0,  for  dynamic  analysis, 

=  1,  for  adjoint  analysis. 

SOLNO  -  an  integer  keeping  a  running  acount  of  the  solution  num¬ 
bers  at  different  time  steps  during  dynamic  analysis. 

DLDB  -  derivatives  of  the  integrands  of  the  cost  and  constraint 
functionals  with  respect  to  b  (design  parameters) . 


DFDB 

DPHDB 

DZIDB 


derivatives  of  the  equations  of  motion,  constraint 
equations  and  spring-damper  relations  with  respect  to  b. 


PB 


-  elements  of  P  matrix  in  equations  of  motion. 
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DNDB  -  derivatives  of  the  initial  solution  variables  with  res- 

,  /3V\ 

pect  to  b,  . 

SUBROUTINE  HALS(NP1,YY, JSIZ) 

The  subroutine  HALS  evaluates  the  integrands  of  the  cost  and  con¬ 
straint  functionals  at  each  time  step. 

SUBROUTINE  INFUNC ( TMAX , ITR , EP I , NFCI , NCV , JS IZ , YY , DLDB , 

DPHDB , DZIDB , DPDFZT , DNDB , DGDB , DGDZ , DPDB , PB ,  ALM) 

The  subroutine  INFUNC  calculates,  through  the  subroutine  AHLPSI, 
the  sensitivity  coefficients  for  the  inequality  functional  constraints 
(Eqs.  (3.10)). 

ITR  -  running  number  of  optimiaztion  iterations 

DGDB  -  derivatives  of  the  non-integral  parts  of  the  cost  and 
constraint  functionals  with  respect  to  b 

DGDZ  -  same  with  respect  to  initial  state  variables 

DPDB  -  derivatives  of  P  matrix  with  respect  to  b. 

SUBROUTINE  TEQFUN(TMAX, ITR,EPE,NFCE ,NFCET, NCV, JSIZ ,YY,DLDB, 

DFDB , DPHDB , DZIDB , DPDFZT , DNDB , DGDB , DGDZ , DPDB , PB , ALM) 

The  subroutine  TEQFUN  calculates  design  sensitivity  coefficients 

for  the  transformed  equality  functional  constraints  (Eqs.  (3.7)). 

SUBROUTINE  DPARMC(TMAX,ITR,EPDI,NDCI ,NCV) 

The  subroutine  DPARMC  calcualtes  design  sensitivity  coefficients 

for  the  design  parameter  constraints  (Eqs.  (3.8)). 

SUBROUTINE  AHLPSI (NCV , ITR , TMAX , ISO , DLDB , DFDB , DPHDB ,DZIDB , 

DPDFZT ,  DNDB ,  DGDB ,  DGDZ ,  DPDB ,  PB ,  YY ,  JS  IZ ,  ALM) 


The  subroutine  AHLPSI  calculates  design  sensitivity  coefficients 


for  the  constraints  and  cost  functionals,  with  the  help  of  sub¬ 


routines  ADJONT,  PBFN,  ADPDB,  and  DGDBZ. 

ISO  -  running  number  of  the  functional  constraint  treated. 

SUBROUTINE  NEWB (EPS ,NCV, STEP, ITR.NFCV) 

The  subroutine  NEWB  computes  the  matrix  of  Eq.  (4.29), 

1  2 

solved  Eqs.  (4.30)  and  (4.31)  for  y  and  y  ,  and  computes  the  nec¬ 
essary  design  changes  6b  given  by  Eqs.  (4.33)  through  (4.35). 

STEP  -  desired  reuction  in  cost  functional 
NFCV  -  number  of  total  functional  constraint  violations. 
SUBROUTINE  DYNANL (IS , IECRO , ITROB , TMIN ,  TMAX,TSTEP ,HMIN, 

HAMX , H , DLDB , DFDB , DPHDB , DZIDB , JSIZ , YY , PB , ALM , LIN , TLIMIT , . 

EPS , IFLAG1 , SOLNO , S0LN01 , ISO , DNDB) 

Subroutine  DYNANL  is  described  in  detail  in  Section  5.3. 

IS  -  a  logic  variable, 

=  0,  for  static  code  generation, 

=  1,  for  static  solution, 

=  2,  for  dynamic  code  generation, 

=  3,  for  dynamic  solution. 

LIN  -  a  user-supplied  flag  that  limits  the  number  of  cor¬ 
rector  iterations  to  1  if  all  equations  are  linear, 

=  0,  for  linear  equations 
=  1,  for  nonlinear  equations. 

S0LN01  -  an  integer  keeping  a  running  account  of  the  solution 

numbers  at  different  time  steps  during  adjoint  analysis. 

SUBROUTINE  ADJONT ( ISO , IFLAG1 , S0LN01 , DLDB , DFDB , DPHDB , DZIDB , 

JSIZ ,YY ,PB , ALM, DNDB) 


94 


The  subroutine  ADJONT  works  in  a  similar  manner  as  ADAMS2  to 
solve  the  adjoint  equations. 

FUNCTION  AHSM(I) 

The  function  subprogram  AHSM  determines  the  maxima  of  the  integrands 
of  the  transformed  equality  constraints,  which  are  used  in  TEQFUN. 

FUNCTION  FUNPSE(I.TMAX) 

The  function  subprogram  FUNPSE  evaluates  the  values  of  the  func¬ 
tional  constraints  and  the  cost  functional. 

SUBROUTINE  PBFN(PB) 

The  subroutine  PBFN  calculates  terms  of  the  P(b)  matrix. 

SUBROUTINE  ADPDB (DPDB , PB) 

The  subroutine  ADPDB  calculates  derivatives  of  the  P(b)  matrix 
with  respect  to  design  parameters. 

SUBROUTINE  DGDBZ (DGDB , DGDZ , ISO) 

The  subroutine  DGDBZ  calculates  the  derivatives  of  the  non-inte¬ 
gral  parts  of  the  cost  and  constraint  functionals. 

SUBROUTINE  DLDFB (  TIMP  U ,  I S  0 ,  YY ,  J  S I Z , DLDB , DFDB , DPHDB , DZIDB ) 

3L  3f  3<f> 

The  subroutine  DLDFB  calcualtes  the  derivatives  ttt  >  rr*  >  XT’  * 

o  b  do  d  b 

3  £ 

and  Tj-g-  for  sensitivity  analysis. 

SUBROUTINE  S01NEW(IS0, JSIZ,HB,PB,AB,ALM,A,T,TIME,ZMAX) 

The  subroutine  S01NEW  calculates  the  right-hand  sides  of  the 
corrector  euqations  for  the  adjoint  system. 

HB  -  time  step  of  dynamic  analysis 

AB  -  transformed  gear  coefficients  of  dynamic  analysis  used 
in  adjoint  analysis 

A  -  current  transformed  gear  coefficients 
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T  -  current  time  step 

TIME  -  current  time 

ZMAX  -  total  time  interval  of  dynamic  analysis. 


SUBROUTINE  ADLDZ ( ISO, JSIZ, TIME, ZMAX) 

The  subroutine  ADLDZ  calculates  the  derivatives  of  the 

integrands  of  the  cost  and  constraint  functionals  with  repsect  to  sol¬ 
ution  variables  of  the  transient  analysis. 

SUBROUTINE  INTERP (Y , N , TA1 , NP3 , TMAX , JSTAT1) 


The  subroutine  INTERP  computes  interpolated  values  of  the  depen¬ 
dent  variable  Y(I,1)  and  its  time  derivatives  and  replaces  the  previous 
values.  The  interpolation  is  to  the  point  TOUT  and  uses  the  Nordsieck 
history  array  Y  as  follows: 


and 


where 


JSTATl 

Y(I,1)  =  l  Y(I,J+1)*S**J, 

J=0 

JSTATl 

Y(I,2)  =  l  J*Y(I,J+1)*S**(J-1), 

J=1 


S  =  (TOUT  +  TA(NP3)  -  TMAX) /HA, 


TA(NP3)  is  ah  old  time  instant  during  backward  integration  and  HA  is 
the  time  step  of  dynamic  analysis. 

SUBROUTINE  VARSET(YY) 

The  subroutine  VARSET  sets  the  variable  YY  to  be  used  in  DIFSUB 
for  numerical  integration 


SUBROUTINE  S03000 (NB ,NJ , IB , JT , NSD , IBSD , NPSOR,NPOSC ,NG ,NT , 
JSIZ , ITF ,NPSS ,NSDT , IBSDT , IMODE) 
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The  subroutine  S03000  sets  up  pointers  of  the  sparse  matrix  non¬ 
zero  position. 

NPOSR  -  pointer  for  the  row  number  of  each  nonzero  entry 

NPOSC  -  pointer  to  the  column  number  of  each  nonzero  entry 

NG  -  an  integer  variable  keeping  a  running  index  of  all 
the  row  and  column  pointers  to  the  nonzero  entries 
in  the  Jacobian  matrix 

NT  -  a  vector  giving  values  of  NG  after  execution  of  diff¬ 
erent  segments  of  S03000 

ITF  -  a  flag  for  keeping  10  nonzero  positions  for  a  revolute 
joint  and  14  for  a  translational  joint 

NPSS  -  a  vector  storing  different  NG's  at  the  beginning  of  each 
block  of  nonzero  positions  for  a  body  in  the  Jacobian. 

SUBROUTINE  S01000 (T , A,H, JSIZ , IMODE ,YY , 1FLAG1 , PB , AB , ISO) 

The  subroutine  S01000  mainly  updates  the  terms  of  the  Jacobian 

matrix  and  the  right  hand-side  terms  of  the  corrector  equation  (see 

Section  5.3  for  further  details). 

SUBROUTINE  INVERT (G , JSIZ , NP2 , IP , JU , JC , IXL , IXH , IPP , ICNT , IPR, 

IPC , IPRI , IVA, IVL , IVU , IC ,IU  ,JA,IRP,IRL, IWSR , ICP , ICL , IWSC , 

IEWM , IRWC , AT , MAXS ,MAXN , IB ,NP0S , IMODE , ITF ,NPSS , IFLAG1 , IPSAV, 

JAS AV , IKNT , KFLAG) 

The  subroutine  INVERT  orders  a  matrix  optimally  and  generates  and 
stores  sparse  matrix  codes  for  LU  factorization.  Provisions  have  been 
made  in  the  program  to  call  it  in  DYNANL  whenever  necessary  during  tran¬ 
sient  and  adjoint  analyses. 

SUBROUTINE  DIFSUB(N,TIME,T,  TMIN,TMAX,EPS,YMAX, ERROR, KFLAG, 

J START ,  MAXDER , G , IVPTR , NERR , NRR , NAL , LIN , CS , DY , TDIV ,NC0NV , 


NPOS , JA , IP , IVA , IC , IU , JC , JU , DI , CC , U , IVL , I VU , CL , IPR ,  IPC , 


IPRI , TLIMIT , IW1 , TNEW , NQ , Y , YS , IFLAG1 ,  PB ,  AB ,  ISO ,  A ,  HA , 

ZMAX,  TA,NP3 , JSTAT1 , PTEST , IKNT) 

The  subroutine  DIFSUB  is  the  numerical  integration  routine  based  on 
the  gear  algorithm  (Chapter  II)  and  is  of  extreme  importance.  For  an 
original  listing  and  explanation  of  the  program  and  variables,  refer¬ 
ence  [7]  is  recommended. 

Explanations  of  the  sparse  matrix  variables  in  subroutines  INVERT 
and  DIFSUB  are  available  in  reference  [61] . 


CHAPTER  VI 


APPLICATIONS  AND  NUMERICAL  RESULTS 

6.1  Introduction 

From  the  structure  of  the  computer  program  DADS  described  in  the 
previous  chapter,  it  is  clear  that  large  classes  of  mechanical  design 
problems  can  be  handled  through  its  implementation.  To  test  the 
program,  two  example  problems  have  been  considered  in  the  following 
sections:  (1)  a  slider-crank  mechanism  (as  described  in  Chapter  II; 

see  Figure  2.5),  and  (2)  a  trip  plowshare  mechanism. 

The  slider-crank  mechanism  is  often  used  in  engines  and  machiner¬ 
ies  [71]  (see  Chapter  II)  and  thus  it  is  quite  familiar.  The  trip 
plowshare  mechanism  treated  is  a  simplified  version  of  the  2500  semi¬ 
integral  spring-reset  plow  which  is  in  production  at  John  Deere 
(see  Section  6.3).  The  importance  of  such  mechanisms  needs  little 
description. 

6.2  The  Slider-Crank  Mechanism 

The  slider-crank  mechanism  to  be  considered  here  is  described  in 
Chapter  II  (Subsection  2. 1.2. 5).  Figure  2.5  shows  the  approximate 
initial  position  of  such  a  mechanism  with  one  spring-damper  pair. 

Link  1  is  ground ,  link  2  is  the  crank  shaft,  link  3  is  the  connecting 
rod  or  coupler,  and  link  4  is  the  piston  or  simply  slider.  A  spring- 


damper  pair  is  attached  between  link  4  and  ground  (Figure  2.5).  There 
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is  a  translational  joint  (type  2)  between  bodies  4  and  1  and  one 
revolute  joint  (type  1)  between  each  of  the  following  pairs  of  bodies: 

1  and  2,  2  and  3,  and  3  and  4.  Gravitational  forces  are  excluded  from 
the  present  simulation  of  the  mechanism. 

The  two-dimensional  ADAMS  program  [14]  has  been  implemented  through 
the  subroutine  ADAMS2  to  obtain  a  static  equilibrium  configuration  and 
to  determine  the  subsequent  transient  motion.  The  results  of  this 
analysis  are  used  later  in  sensitivity  analysis  and  optimization. 

A  symbolic  listing  of  the  nonzero  positions  of  the  Jacobian  matrix 
for  this  example  problem  is  given  in  Fig.  2.7.  An  explanation  of  the 
nonzero  positions  can  be  obtained  in  reference  [14]. 

6.2.1  Formulation  of  the  Optimal  Design  Problem 
By  virtue  of  its  movement,  a  radial  slider-crank  mechanism  exerts 
a  force  on  ground  through  the  crank-bearing  and  the  wrist-pin  guide 
(such  as  cylinder  wall  in  an  automotive  piston-type  engine) .  It  is 
desirable  to  keep  these  "shaking  forces"  within  bounds.  It  is  also 
desirable  to  put  an  upper  bound  on  the  angular  velocity  of  the  crank 
at  the  final  instant  T  of  the  time-interval  [0,T]  under  consideration. 
The  cost  function  is  chosen  to  be  twice  the  maximum  energy  stored  in 
the  spring  during  the  interval  of  motion.  The  following  design  param¬ 
eters  are  considered  (see  Chapter  II) : 

b^  =  The  spring  constant  of  the  spring 

b^  =  Height  of  the  points  of  attachment  of  the  spring 

b^  =  Half  of  the  length  of  the  uniform  coupler. 
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With  the  notations  of  Chapters  II  and  III,  the  optimal  design 
problem  is  stated  as  follows:  minimize 


J  = 


max  b,  (& ,  -  £  n  ) 
0<t  <T 


(6.1) 


subject  to  the  equations  of  motion  (2.33),  the  equations  of  con¬ 
straints  on  motion  (2.35),  the  spring-damper  relations  (2.36),  the 
initial  conditions  of  the  form  of  Eqs.  (3.4),  and  (3.5),  the  func¬ 
tional  constraints 

|  y  2 1  £  P2(max)  *  Oft  <  T  (6.2) 

<£2(T)  _<  (jl2 (max)  (6.3) 


and  the  design  parameter  constraints 

bilbi<bJ  >  i-1,2,3  (6.4) 

In  Eq.  (6.2)  p„  is  the  y-component  of  reaction  forces  acting 
'  3$y  3$2 

on  body  2  at  joint  1,  because  from  Eq.  (2.19)  r-r—  =  r--—  =  -  1, 

a$2  2  9Y2 

and  so  -  y2  =  P2  in  Eq’  ^2*  2^)  • 

For  the  sake  of  simplicity,  only  the  constraints  on  the  vert¬ 
ical  component  of  shaking  forces  at  the  crank  bearing  have  been 
considered  here. 

After  the  introduction  of  an  artificial  design  parameter  b^ 
[24,35,37]  and  the  integral  functional  forms  in  the  conventional 
way  [24],  the  problem  can  be  reformulated  as:  minimize 


*■ 


*0  =  b4 


(6.5) 
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subject  to  Eqs.  (2.33),.  (2.35),  (2.36),  (3.4)  and  (3.5),  the 
constraints 

<T  y2  ~y2(max)>dt  =  0  (6*6) 

<y2-y2(max)>t  =  0  (6'7) 

fT  2 

^ 3  =  ]o  <b1(£1-£1Q)Z  -  b4>dt  =  0  (6.8) 

*4  £  *2^  “  ^2 (max)  -  0  (6’9^ 


and  the  design  parameter  constraints  of  Eq.  (6.4). 

The  symbol  <q  (t£>  used  above  has  the  following  meaning  (see 
Chapter  III) : 


<n  (t)> 


for  n  (t)  >_  0 
for  n  (t)  <  0 


(6.10) 


This  formulation  corresponds  to  the  general  formulation  of  the 
optimal  design  problem  stated  in  Chapter  III.  After  specification 
of  initial  numerical  data,  the  DADS  program  described  in  Chapter  V 
can  be  implemented  to  obtain  a  solution. 

6.2.2  Sensitivity  Analysis 

In  problems  of  optimal  design,  sensitivity  analysis  constitutes 
a  principal  part  of  the  work  to  be  done.  When  sensitivity  coeffic¬ 
ients  are  obtained,  an  iterative  optimization  algorithm  can  be 
applied  to  determine  an  optimal  solution.  In  the  present  class  of 
problems,  subroutines  RELATE,  DNUDB,  DGDBZ ,  ADLDZ,  and  DLDFB  are  the 
major  user-supplied  subprograms  for  sensitivity  analysis.  Among 
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them  DGDBZ  and  ADLDZ  are  often  relatively  simpler,  because  they 

l 

depend  solely  on  the  forms  of  the  functional  constraints.  The  rela¬ 
tions  entering  into  the  subroutines  RELATE,  DGDBZ,  ADLDZ,  and  DLDFB 
for  the  present  problem  (with  1  second  interval)  and  the  plow-share 
mechanism  are  given  in  the  Appendix.  Subroutines  DGDBZ  and  DLDFB 
are  used  in  AHLPSI  to  calculate  sensitivity  coefficients. 


6.2.3  Numerical  Results 

Initial  estimates  of  the  parameters  for  the  slider-crank 
mechanism  under  consideration  are  given  in  Table  6.1.  The  units 
used  are  inch,  pound-force,  and  second.  Two  time  intervals  are 
treated:  [0.0, 1.0]  and  [0.0,  2.0].  During  transient  analysis  a 

constant  counter-clockwise  torque  of  100  in/lbf  is  applied  to  link 
2.  Fig.  6.2  gives  the  tableau  of  nonzero  positions  of  the  Jacobian 
matrix  for  this  problem. 

To  examine  the  correctness  of  Eqs.  (3.40)  and  (3.41)  for 

ij)  g 

sensitivity  analysis,  the  sensitivity  coefficeints  l  ,  6e{l,2,3,4}, 

*0 

and  *•  are  first  calculated  for  the  2  second  interval,  with  initial 
I  T 

estimates  b  =  [1.0,0.5,10.0,14.5]  for  the  design  parameters,  for 

L  T  U  T 

lower  bounds  b  =  [0.8 ,0.2,5 .0]  and  upper  bounds  b  =  [1.5,0.8,16.0] 

on  the  first  three  design  parameters,  and  with  ^(max)  =  10.0  and 

^2 (max)  =  The  cost  and  constraint  functionals  are  evaluated 

for  0.1%  and  1.0%  perturbations  of  the  first,  second,  and  fourth  de¬ 


sign  parameters  (in  conjunction  with  0.01%  and  0.1%  changes  in  the 

third)  and  the  corresponding  changes  in  the  functionals  are  exam- 

*0 

Table  6.2  gives  the  sensitivity  coefficients  l  and  l  p 


ined . 


Table  6.1 


Initial  Estimates  of  the  Parameters  for  the  Slider-Crank 
Mechanism  (Inch,  Pound-Force,  Second). 


LINK  DESCRIPTION 


M(l) 

=  1. 

M(2)  =  4. 

M(3)  =1. 

M(4) 

=  3. 

JIN(l) 

=  1. 

JIN(2)  =10. 

JIN(3)  =4. 

JIN(4) 

=  2. 

X(l) 

=  0. 

X(2)  =1. 

X(3)  =8.667 

X(4) 

=  17.32 

Y(l) 

=  0. 

Y(2)  =0. 

Y(3)  =5.25 

Y(4) 

=  .5 

PHI(l) 

=  0. 

PHI  (2)  =1.5708 

PHI(3)  =-.5236 

PHI (4) 

=  0. 

FX(1) 

=  0. 

FX(2)  =0. 

FX(3)  =0. 

FX(4) 

=  0. 

FY(1) 

=  0. 

FY  (2)  =0. 

FY  (3)  =0. 

FY(4) 

=  0. 

TQ(1) 

=  0. 

TQ(2)  =0. 

TQ(3)  =0. 

TQ(4) 

=  0. 

JOINT  DESCRIPTION 

IB(1, 1) 

=  1 

IB(1,2)=  2 

IB(1,3)=  3 

IB(1 ,4) 

=  4 

Xl(l) 

=  0. 

Xl(2)  =9.0 

XI  (3)  =10. 

XI  (4) 

=  0. 

Yl(l) 

=  0. 

Yl(2)  =0. 

Y1  (3)  =0. 

Yl(4) 

=  0.5 

IB (2,1) 

=  2 

IB(2,2)=  3 

IB(2,3)=  4 

IB (2, 4) 

=  1 

X2(l) 

=  -1. 

X2(2)  =-10. 

X2(3)  =0. 

X2  (4) 

=  0. 

Y2(l) 

=  0. 

Y2(2)  =0. 

Y2  (3)  =0. 

Y2(4) 

=  1.0 

SPRING-DAMPER  DESCRIPTION 

IBSD(1 ,1) 

=  1 

XF1(1)  =35. 

YF1(1)  =0.5 

IBSD(2,1) 

=  4 

XF2(1) 

=  2. 

YF2(1)  =0. 

SK(1)  =1. 

DC(1) 

=  1. 

SDL(l) 

=  15. 

68  SDL0(1)=  15.68 

SDF(l)  =0. 

Table  6.2 


Sensitivity  Analysis  Results  for  the 


105 


% 


4J 

3 

§* 

o 

u 


I 

cu 

S 

a 

rt 


QJ 

3 

rH 

CO 


CM 

e- 

a 

3  I 

Y 

I  0 

I 

\  N 

n 

H 

rH 

3- 

a 

3  1 

V  1  0 

IAN 

n 

m 

-1 

■ 

•o 

o 

na 

9 

rH 

o 

1 

X 

rH 

EH 

H 

in 

m 

<3* 

m 

*0 

IE 

o 

• 

o 

o 

rH 

o 

o 

• 

o 

rH 

» 

« 

■ 

O 

o 

o 

o 

o 

o 

o 

o 

■ 

-1 

I  : 

'g 

o 

■ 

rH 

• 

tH 

V-' 

*»n 

BI 

m 

m 

**3* 

il 

HI 

o 

o 

o 

rH 

o 

o 

o 

rH 

Q 

o 

o 

o 

o 

o 

o 

o 

o 

,nY 

•  • 

in 

CO 

1 

m 

u 

Q) 

•  ►a 

m 

m 

m 

m 

rH 

m 

m 

m 

m 

vO 

4J 

• 

• 

• 

• 

• 

• 

• 

• 

• 

0) 

*<r 

<r 

*■3* 

s 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

nj 

L  . 

u 

■M 

-1 

CO 

fh 

*o 

rH 

A 

60 

■ 

X 

•H 

<r 

m 

m 

CO 

aRH 

<r 

a) 

<o 

B 

O 

o 

o 

rH 

o 

o 

O 

rH 

Q 

<U 

x: 

M 

o 

o 

o 

o 

o 

o 

O 

o 

4J 

CM 

H 

IM 

o 

O 

o 

rH 

rH 

Q) 

cn 

rO 

X 

X 

■U 

o 

o 

o 

rH 

o 

o 

o 

rH 

o 

• 

• 

% 

• 

• 

• 

• 

• 

• 

a 

i  ° 

o 

o 

o 

o 

o 

o 

o 

o 

U 

■■ 

CO 

w 

1 

CO 

N 

1 

■ 

o 

o 

rH 

rH 

cO 

CM 

rO 

pi 

X 

■U 

<o 

E3 

o 

m 

o 

o 

o 

in 

o 

o 

•H 

C 

i  o 

o 

o 

o 

o 

o 

o 

o 

o 

M 

. 

CM 

v 

o 

o 

rH 

rH 

rH 

‘O 

x 

x 

o 

rH 

o 

o 

o 

rH 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

Table  6.3  (Cont’d) 


6 


107 


Be  {3,4},  and  Table  6.  3  gives  other  computational  results  in 
compact  form.  The  first  and  second  constraints,  however,  remain 
unviolated . 


In  these  computations  the  integrands  of  t(i ^ »  and  ^3  have 

2  2  2 

been  normalized  by  dividing  by  ^(max)’  ^2 (max)’  and  ^4’  resPect_ 
ively.  It  is  observed  from  the  Table  that  terms  of  the  form 
Jl  6b  match  satisfactorily  with  the  corresponding  Alp's.  The  slight 
discrepancies  are  attributed  to  various  approximations  and  lineariz¬ 
ations  at  many  steps  of  both  transient  and  sensitivity  analysis  and 
relatively  coarse  time  grids  for  numerical  integrations  and  trans¬ 
formation  of  pointwise  constraints  to  functional  forms.  However, 
they  need  not  affect  the  optimization  process. 

In  carrying  out  the  optimal  design  algorithm,  a  design  reduct¬ 
ion  ratio  of  3%  is  used  to  compute  the  stepsize  in  the  first  iter- 

L  U  T 

ation.  The  bounds  b  and  b  are  taken  to  be  [0.8 ,0.2, 9.0]  and 

T  I 

[1.5,0.8,12.0]  ,  respectively.  The  initial  estimate  b  (including 


artificial  design  parameter)  is  taken  as  [0.8244,  0.5,  9.0,  11.12] 

and  p 2 (max)  and  ^ 2 (max)  are  kept  unchanged • 

At  the  starting  design,  | | 6b  ||  =  0.1248  and  | | 5b  ||  =  0.2307, 

constraints  3  and  4  are  violated,  and  the  constraint  on  the  design 

parameter  b^  is  tight.  After  6  iterations,  as  the  algorithm 

approaches  the  optimum  design,  |  |  6b'*' |  |  reduces  to  0.8404x  10 
2  “1 

|  1 6b  ||  reduces  to  0.2669x  10  ,  and  the  constraints  on  b^,  b^,  and 

b^  become  tight.  The  final  optimum  results  are  given  in  Table 


6.4. 


Table  6.4 


Optimum  Results  for  the  Slider -Crank  Mechanism 
for  the  Time  Interval  of  2  Seconds 


Real  cost  function  J 
Lower  bounds  on  b 
[Upper  bounds  on  b 


max  b1  (£.  - 
0 _<t<2 

[0 .8 , 0. 2 , 9 .  Oj 
[1.5,0. 8,12. 0]T 


Starting  Values 


Optimum  Values 


8.2440  x  10 


-1 


8.0000  x  10 


-1 


5.0000  x  10 


-1 


5.0000  x  10 


-1 


9.0000 


9.0000 


1.1159  x  10 


1.00115  x  10 


Next,  the  design  problem  is  considered  for  the  time  interval 


of  1  second  with  a  different  set  of  data.  For  this  case  a  de¬ 
sign  reduction  ratio  of  10%  is  used  in  computing  the  step  size 
in  the  first  iteration.  The  bounds  b1  and  b^  are  taken  to  be 
[0.8,  0.2,  5.0]1  and  [1.5,  0.8,  16. 0]1,  respectively.  The  init- 


I  T 

ial  estimate  b  is  taken  as  [1.0,  0.5,  10.0,  1.0]  and  p 


and  $2 (max)  are  ta^en  as  5*0  and  0.3,  respectively. 


At  the  starting  design, 


and  constraint  3  is  violated. 


fib  I  I  =  0.0,  I j fib  . . 

In  the  second  iteration 


2  (max) 


=  0.1106x10 
fib1 1  I 


2.  —I 

0.7348,  |  |  fib  ||  =  0.3814x10  ,  and  again  constraint  3  is 


violated.  After  18  iterations,  as  the  algorithm  approaches  the 
optimum  design,  |  |  fib"1 1  |  =  5.567x  10  ^,  ||fib^||  =3.018x10^, 
and  (lower)  constraints  on  b.^  and  b2  remain  tight.  The  optimum 
results  are  given  in  Table  6.5. 


6.3  The  Plow-Share  Mechanism 

The  transient  dynamic  response  and  design  sensitivity  of  a 
spring-rest  plow-share  mechanism  is  determined  to  illustrate  the 
flexibility  of  the  DADS  computer  program  .for  systems  with  inter¬ 
mittent  motion.  Fig.  6.1  shows  the  approximate  initial  position 
of  such  a  mechanism.  For  a  system  of  this  nature  any  type  of 
closed  form  solution  is  beyond  consideration. 

The  model  consists  of  6  moveable  rigid  bodies,  identified  as 
follows:  body  2  -  plow-share  and  standard;  body  3  -  lower  link; 


Joint  3 


Approximate  Initial  Configuration  of  the  Plow-Share 


Table  6.5 


Optimum  Results  for  the  Slider  Crank 
Mechanism  for  the  Time  Interval  of  1  Second 


-  12 
Real  Cost  Function  J  =  max  b.  (2..  -  SL  ) 

0<t<l  0 


Lower 

Bounds 

on  b  =  [0.8, 

0.2, 

5 . 0] T 

Upper 

Bound  s 

on  b  =  [1.5, 

0.8, 

16.0]' 

Starting 

Values 

Optimum 

Values 

h 

1.0000 

8. 0000 x  10"1 

b2 

5 . 0000  x  10"1 

5. 0000 x  10_1 

b3 

1 . 0000 x  101 

9.7566 
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body  4  -  rear  toggle  link;  body  5  -  front  toggle  link;  body  6  - 
U-bolt;  and  body  7  -  combined  plow  frame  and  tractor.  Body  1 
(ground)  is  rigidly  connnected  to  the  inertial  reference  frame. 
Various  pinned  connections  (revolute  joints)  are  shown  in  Fig. 

6.1.  The  entire  system  is  moving  to  the  right  at  approximately 
2  meters  per  second,  modeled  by  a  horizontal  translational  joint 
between  ground  (body  1)  and  the  combined  plow  frame  and  tractor 
(body  7).  A  linear  spring-damper  combination  is  connected  bet¬ 
ween  the  U-bolt  and  rear  toggle  link.  In  addition,  5  contact 
points,  identified  by  the  letters  A-E,  represent:  A  -  contact 
between  the  U-bolt  and  main  frame;  B  -  contact  between  the  shank 
and  lower  link;  C  -  contact  between  the  lower  link  and  main 
frame;  D  -  contact  between  the  front  and  rear  toggle  links;  and 
E  -  contact  between  the  plow-share  tip  and  rock  embedded  in  the 
ground  (body  1) . 

Contacts  are  simulated  by  attaching  modified  linear  spring- 
damper  combinations  between  contacting  bodies.  The  modification 
consists  of  setting  the  spring  and  damping  coefficients  to  zero 
in  a  continuous  manner  as  parts  break  contact  and  to  their  max¬ 
imum  values  as  contact  is  made.  This  is  accomplished  by  intro¬ 
ducing  spring  and  damping  coefficients  as  state  variables  and  mak¬ 
ing  them  functions  of  spring  length.  A  three-parameter  model  is 
chosen  because  of  the  wide  range  of  spring  and  damper  character¬ 
istics  it  is  capable  of  representing.  The  equations  are  of  the 
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form 

b.  . 

k.  .  =  xTfd  -  exp  [-( <  ±  (a  -  £  n)  >  /e. .)  13  ]) 

1J  1J  O'  lj 

b.  . 

Cij  =  Cij(1  '  EXP[-(<  1  Ci -*<,)>  /0ij)  13 1)  (6.11) 

where  and  are  the  maximum  values  of  the  spring  and  damping 
coefficients  for  the  pair  connecting  i^  and  j bodies.  The  ex¬ 
pression  <  ±(£-£q)>  equals  zero  if  ±  (£  -£^)  0;  otherwise  it 

equals  ±  (£  ~  £  q) •  The  *  +  *  or  '  -  '  sign  is  selected  depending 
upon  whether  the  coefficients  are  to  increase  or  decrease  as  the 
spring  length  increases.  The  parameter  0  determines  the  inter¬ 
val  of  spring  travel  over  which  the  coefficients  change  and  b^ 
determines  the  shape  of  the  curve.  Fig.  6.2  depicts  various 
shapes  for  combinations  of  b^  and  9  „ .  For  the  5  contact  points 
(stops)  in  the  present  problem  the  1  -  '  sign  is  chosen  in  equation 
(6.11). 

It  should  be  observed  here  that  by  varying  the  free  spring 

length  (undeformed),  K^,  ^21’  ^21’  anc*  ®21  ^or  t*ie  s*xth  spring, 
various  models  of  the  underground  rock  can  be  obtained .  In  the  pre¬ 
sent  analysis  the  frictional  forces  acting  between  the  plow  tip  and 
the  rock  have  been  taken  as  a  continuous  function  of  the  spring 
forces  and  have  been  fed  into  the  program  through  subroutine 
FOREXT  (see  Chapter  V)  after  transforming  them  to  the  equivalent 
system  at  the  center  of  mass  of  the  second  body. 

Owing  to  the  presence  of  the  nonlinear  spring -dampers,  the 
Jacobian  matrix  changes  very  rapidly  in  a  neighborhood  of  certain 


MAGNITUDE  1  . .  MAGNITUDE 


Figure  6.2  Dependence  of  Functions  on  B  and  Theta. 
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time  instants.  Provisions  have  been  made  in  the  program  to  regen¬ 
erate  the'  sparse  matrix  codes  at  those  time  instants,  by  introduc¬ 
ing  several  error  tests  and  flags  in  subroutine  DIFSUB  (see  Chapter 
V) .  The  error  tests  measure  the  ratios  of  maximum  valued  elements 
of  Jacobian  to  the  pivots.  When  a  pivot  tends  to  zero,  the  ratio 
becomes  very  large  and  a  fresh  code  generation  is  demanded.  For 
the  estimation  of  the  bound  and  related  discussions,  references 
[73,74]  are  recommended. 

6.3.1  Numerical  Results  (Dynamic  Analysis) 

Initial  estimates  of  the  parameters  for  the  spring  reset  plow¬ 
share  mechanism  are  given  in  Table  6.6.  The  units  used  are  meter, 
kilogram,  second,  and  Newton.  For  this  problem,  the  values  of  the 

parameters  b..  and  9..  are  the  same  for  all  five  nonlinear  springs. 

IJ  -*-J 

They  are  1.2  and  8.617739x  10  respectively.  The  corresponding 

-5 

spring  travel  is  approximately  1.538x  10  to  attain  99%  of  the 

maximum  values  of  the  spring-damper  coefficients.  The  maximum 

values  of  the  coefficients  of  the  5  nonlinear  spring-damper  pairs 

are  (l.OxlO6,  6.2xl04),  (l.OxlO6,  6.0xl04),  (l.OxlO6,  6.2xl05), 
6  4  5 

(1.0x10  ,  6.2x10  ),  and  (1. Ox  10  ,  0.0),  respectively.  These 
are  entered  into  the  program  through  subroutine  INDATA  (see  Chapter 
V) .  During  transient  analysis,  the  sparse  matrix  codes  are  regen¬ 
erated  ten  times. 

Initially  the  combined  plow  and  tractor  system  shown  in  Fig. 

6.1  is  moving  to  the  right  at  2  meters  per  second,  along  a 


Initial  Estimates  of  the  Parameters  for  the  Spring-Reset 
Plow-Share  Mechanism  (Meter,  Kilogram,  Second,  Newton) 
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Table  6.6  (cont'd) 
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horizontal  translational  joint.  The  tip  of  the  plow  makes  con¬ 
tact  with  the  rock  at  time  =  0.0  seconds.  As  shown  in  Fig.  6.3, 
the  tip  breaks  contact  with  the  rock  at  about  0.1  seconds,  but 
fails  to  clear  it  and  comes  to  contact  again  between  0.22  and 
0.32  seconds.  The  contact  force  imparts  an  angular  velocity  to 
the  plow-share,  causing  it  to  move  rearward  and  upward  (see  Figs. 
6.4  and  6.5).  This  motion  drives  the  toggle  links  upward, 
bringing  spring  1  into  tension  (see  Fig.  6.6).  The  U-bolt  and 
lower  link  come  into  contact  with  the  plow  frame  (contact  points 
A  and  C)  at  0.11795  and  0.39581  seconds,  respectively.  Contact 
at  B  between  the  standard  and  lower  link  (stop  2)  is  broken  at 
0.32384  seconds  and  this  event  coincides  with  the  loss  of  con¬ 
tact  of  the  tip  with  the  rock. 

Contact  at  C  between  the  lower  link  and  frame  stops  upward 
movement  of  the  plow-share  and  the  reset  cycle  begins.  Stored 
energy  in  the  spring  rapidly  collapses  the  toggle  links.  This 
action,  shown  in  Figs.  6.4  and  6.5,  causes  a  rapid  change  in 
angular  displacement  of  the  plow-share,  with  only  a  small  effect 
on  its  vertical  displacement.  At  0.53509  seconds,  the  toggle 
links  have  reset  (contact  at  D) . 

The  lower  link  and  U-bolt  break  contact  with  the  frame  at 
0.55585  and  0.68854  seconds,  respectively.  It  is  interesting  to 
note  that  the  toggle  action  results  in  the  plow-share  being  brought 
to  within  20°  of  horizontal,  while  its  center  of  mass  is  still 
0.75  meters  above  ground.  The  plow-share  therefore  re-enters  the 


Figure  6.3  Stop  5  Reaction  Between  Plow-Share  Tip  and  Rock. 
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Figure  6.4  Vertical  Displacement  of  Plow-Share  CM. 
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ground  at  a  shallow  angle,  preventing  the  mechanism  from  being 
tripped  again.  Finally,  at  about  0.83174  seconds,  contact  occurs 
at  stop  2  and  the  mechanism  regains  its  approximate  initial  con¬ 
figuration.  Fig.  6.7  gives  the  stop  2  reactions  during  the  en¬ 
tire  time  interval. 

Fig.  6.8  gives  VERSATEC  plots  (snap-shot  pictures)  of  the 
mechanism  at  selected  instants  of  time.  The  whole  process  of  the 
computer  simulation  of  this  dynamic  analysis  takes  about  3  minutes 
30  seconds  of  CPU  time  with  an  IBM-370-168  computer  system. 

6.3.2  Formulation  of  a  Trip-Plow  Optimal  Design  Problem 

It  was  noted  in  the  last  section  that  while  re-entering  the 
ground  the  plow-share  mechanism  should  be  prevented  from  being 
tripped  again.  This  can  be  achieved  if  a  constraint  is  imposed 
such  that  the  plow-share  is  brought  to  horizontal,  to  within  a 
certain  small  tolerance  angle,  during  the  final  phase  of  motion. 

The  cost  function  to  be  minimized  in  the  design  process  is  chosen 
to  be  twice  the  maximum  energy  stored  in  the  spring  during  the  int¬ 
erval  of  motion,  which  is  taken  to  be  0.832  seconds.  The  follow¬ 
ing  design  parameters  are  considered  here: 

b^  =  Diameter  of  the  coil  wire  of  the  reset  spring; 
b 2  =  Diameter  of  the  coil  of  the  reset  spring; 
b^  =  Number  of  turns  in  the  coil. 

The  relation  between  the  undeformed  length  £q  of  the  reset 
spring  and  design  parameters  is  taken  as 
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l 


1 

0 


^1^3  +  0*56604  x  10 


(6.12) 


It  should  be  observed  here  that  since  the  undeformed  length 
of  the  reset  spring  depends  on  b^  and  is  automatically  a 

design  parameter. 

With  the  notations  of  Chapters  II  and  III,  the  optimal  design 
problem  is  stated  as  follows:  Minimize 


J  5 


max 

0<t<0.832 


kV^jiJ)2 


where  [72] 


Gb 


K  =  K 


46 


Sb2b3 


G  being  shear  modulus, 


(6.13) 


(6.14) 


subject  to  the  equations  of  motion,  Eq.  (2.33),  the  equations  of 
constraints  on  motion,  Eq.  (2.35),  the  spring-damper  relations, 

Eqs .  (2.36)  and  (6.11),  initial  conditions  of  the  form  of  Eqs.  (3.4) 
and  (3.5),  the  functional  constraint 


-  <fr2  -  0.174533  <0  ,  0.75  <  t  <_  0.832  (6.15) 

and  the  design  parameter  constraints 

b^1  <  b.  <  bV  ,  i=l ,2,3  (6.16) 

1  -  1  -  1 

For  the  functional  constraint  Eq.  (6.15),  the  maximum  allowable 
inclination  of  the  plow-share  during  the  interval  0.75  _<  t  <_  0.832 
of  re-entering  phase  has  been  taken  to  be  10°. 

As  in  the  case  of  the  slider-crank  mechanism,  with  similar 


* 
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notations,  the  problem  can  be  reformulated  as:  Minimize 


4<0  -  b4  (6.17) 

subject  to  Eqs.  (2.33),  (2.35),  (2.36),  (6.11),  (3.4),  and  (3.5), 
the  constraints 


■■■£ 

'.‘l 


0.832 

<  <  -  d>0  -  0.174533  >  >  dt  =  0 


0.832  Qb 


37"  '  b4  >  dt 


=  o 


(6.18) 

(6.19) 


8b2b3 


and  the  design  parameter  constraint  of  Eq.  (6.16). 

In  Eq.  (6.18)  the  symbol  «  H >>  has  the  following  meaning: 


<  <H>>  = 


0,  for  0<_t<_t-^  =  0.75  and  for  H  <_  0, 


H  ,  for  H  >0 


(6.20) 


Since  the  transient  analysis  itself  is  extremely  complex  for  this 
mechanism,  only  one  regular  functional  constraint  has  been  consid¬ 
ered.  However,  additional  constraints  can  be  treated. 

6.3.3  Modifications  in  Sensitivity  Analysis  Due  to 
Non  Standard  Elements 

Owing  to  the  inclusion  of  nonstandard  equations,  Eq.  (6.11), 
in  the  set  of  spring-damper  relations,  some  modifications  are  re¬ 
quired  in  the  sensitivity  analysis.  The  calculations  and  program¬ 
ming  for  the  additional  elements  in  NPOSR,  NPOSC,  and  G  vectors 
and  the  right-hand  sides  are  routinely  [7,14]  done  through  the 
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subroutines  USET,  US0LV1,  US0LV2,  and  US0LV3  together  with  another 
set  of  subroutines  LSPTRS,  LSSLVl ,  LSSLV2,  and  LSSLV3  (see  Chapter 
V) .  The  solution  vector  of  of  Eq.  (6.11)  is  included 

in  the  extended  vector  and  the  corresponding  I  matrix  of  Eq. 
(2.38)  becomes 


I 


0  0 
1  0 

0 

0 

0 

1 


j(4s  x  4s) 
0 

0 


(6.21) 


0 


(4s  +  NSDV) 
x  (4s  +  NSDV) 


where  NSDV  is  twice  the  number  of  nonlinear  spring-damper  pairs. 

For  the  sensitivity  analysis,  necessary  modifications  are  per¬ 
formed  in  subroutine  S01NEW  for  calculation  of  the  right-hand  sides 

/ 

of  the  adjoint  equations  and  additional  calculations  for  the  values 
3  £ 

of  yj-  are  done  in  DLDFB.  Necessary  adjustements  for  the  dimensions 
of  the  variable  DZIDB  are  also  made.  It  should  be  remarked  here 


♦ 
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that  DLDFB  is  called  after  the  call  of  INTERP,  so  that  the  values 
in  DLDFB  are  calculated  with  the  interpolated  values  of  the  sol¬ 
ution  variables. 

6.3.4  Numerical  Results  (Adjoint 
Analysis  and  Optimization) 

In  carrying  out  the  optimal  design  algorithm,  a  design  re¬ 
duction  ratio  of  1%  is  used  to  compute  the  step-size  in  the  first 

L  U  -2 

iteration.  The  bounds  b  and  b  are  taken  to  be  [0.2x10  , 

0.1x10  0.8x10  ]  and  [0.1x10  0.1,  0.14  x  10^]  ^ ,  respect¬ 

ively.  The  initial  estimate  b'*'  (including  artificial  design  para¬ 
meter)  is  taken  as  [0.56604x  10  0.181361x10  \  0.1x10^, 

4  T 

0.167471x10  ]  .  The  shear  modulus  G  has  been  taken  to  be 

1.86xl01(^  N/m2.  The  static  analysis  and  static  sensitivity  anal¬ 

ysis  results  are  given  in  Table  6.7  and  Table  6.8. 

It  is  interesting  to  note  that  in  the  first  iteration,  all 

1  "“3 

convergence  criteria  are  satisfied  with  |  |  6b  ||  =  0.6209x  10 
2  —3 

and  |  1 6b  ||  =  0.2132x10  and  the  final  optimum  results  are  given 
in  Table  6.9. 


In  order  to  treat  the  problem  in  a  different  way,  the  first 
functional  constraint  is  redefined  such  that  the  tip  of  the  plow¬ 
share  remains  at  least  2"  =  0.0508  m.  above  the  ground  level  dur¬ 
ing  a  certain  interval  in  the  final  phase  of  motion.  This  time 
interval  is  taken  to  be  [0.6,  0.832].  The  corresponding  functional 


'0 


(6.22) 


Static  Sensitivity  Analysis  Results 
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See  Chapter  V  for  the  variable  names. 
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where 

H  =  0.2508  -  (Y(2) +XF1(6)*DSIN(PHI(2)) 

+  YF1.(6)*DC0S  (PHI(2))]  (6.23) 

and 

<<H  » 

The  computer  variables  used  in  Eq.  (6.23)  are  explained  in  Chapter 
V. 

For  this  case  a  design  reduction  ratio  of  5%  is  used  and  the 
bounds  b^  and  b^  are  kept  unchanged.  The  initial  estimate  B*  is 
taken  as  [0.566xl0-2,  0.181361 x  10_1,  O.lxlO2,  0.167471x  104]T. 

In  this  case,  optimum  results  are  obtained  in  the  first  iteration, 
with  ||Sb'*'j|  =  0.0  and  ||Sb2|)  =  0.2011x10  4.  The  final  results 
are  given  in  Table  6.10. 

It  should  be  remarked  here  that  although  the  results  in  Tables 
6.9  and  6.10  are  quite  similar,  the  results  of  Table  6.9  are  prefer¬ 
able,  because  the  constraint  on  inclination  of  the  plow-share  is 
more  important  than  the  one  on  the  height  of  the  plow-share  tip. 

In  order  to  deal  with  a  more  meaningful  problem  of  iterative 
optimal  design,  the  maximum  inclination  of  the  plow-share  during 
the  interval  0.6  _<  t  _<  0.832  of  the  re-entering  phase  is  reduced  to 
0.017453  radians.  Then  Eqs.  (6.15)  and  (6.18)  reduce  to 

-  <f)2  -  0.017453  _<  0  ,  0.6  ^  t  <  0.832  (6.25) 


jo  ,  for  0  <_  t  _<  t^  =  0.6  and  for  H  _<  0 
1h2,  for  H  >  0 


(6.24) 


and 


Table  6.9 


Table  6.10 


Optimum  Results  for  the  Spring-Reset 
Plow-Share  Mechanism  with 
Modified  Functional  Constraint 


Ght  12 

Real  Cost  Function:  J  =  max  - - -  (£  -Sin) 

0  <t_<0 . 8 3 2  8  b,  b  1  U 

-2  -1  IT 

Lower  bounds  on  b  =  [0.2x10  ,  0.1x10  ,  0.8x10  ] 

Upper  bounds  on  b  =  [0.1x10  \  0.1,  0.14x  lO2]'*' 

Starting  Values 

Optimum  Values 

bl 

0.566 x  10~2 

0.567958 x  10"2 

b2 

0.181361 x  10-1 

0.181315 x  10-1 

b3 

0.1 x 102 

0.1  x  102 

b4 

0.167471 x 104 

0.167471 x  104 

(app . ) 
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•  rO .832 

ij>  =  <<  -  <f>  -  0.017453  »  dt  =  0  (6.26) 

J0 

respectively.  All  other  equations  in  the  set  of  Eqs.  (6.12)  to 
(6.20)  remain  unaltered. 

For  this  case  a  design  reduction  ratio  of  0.25%  is  used  and 
the  bounds  and  b^  are  kept  unchanged.  The  initial  estimate 
b1  is  taken  as  [0.56604 x 10-2 ,  0. 181361 x  10_1,  O.lxlO2, 

4  x 

0.16000x  10  ]  .  After  three  iterations  the  convergence  criteria 
are  satisfied  and  the  optimum  results  are  reached.  Tables  6.11 
and  6.12  present  the  pertinent  numerical  results  of  interest  for 


various  iterations. 


Numerical  Results  for  the  Spring-Reset  Plow-Share 
.anism  (with  0.017453  Radians  as  the  Maximum  Allowable 


Table  6.12 
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CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  research,  a  systematic  and  unified  theory  for  obtaining 
a  corresponding  computer  program  for  dynamic  analysis,  static  sensiti¬ 
vity  analysis,  general  design  sensitivity  analysis,  and  optimal  design 
of  large  scale  constrained  dynamic  systems  has  been  developed  and  demon¬ 
strated.  Constrained  mechanical  systems  with  continuous  motion,  as  in 
the  case  of  the  slider-crank  mechanism  of  Chapter  VI,  and  those  with 
intermittent  motion,  as  in  the  case  of  the  plow-share  mechanism  dis¬ 
cussed  in  Chapter  VI,  are  handled  with  ease.  The  stiff  integration 
(Gear)  algorithm  [6,2]  and  sparse  matrix  techniques  [6,10,60,61,62]  have 
been  successfully  employed  in  the  DADS  program  that  is  presented  here. 
The  DADS  program  is  presently  configured  for  two-dimensional  systems, 
but  the  extension  to  three-dimensional  systems  is  theoretically 
straightforward. 

It  was  noted  in  Chapter  V  that  the  DADS  computer  program  is 
capable  of  dealing  with  mechanical  systems  with  discrete  rigid 
bodies  only.  If  necessary,  some  amount  of  def ormability  in  special 
elements  of  mechanisms  can  be  introduced,  with  the  help  of  additional 
spring-damper  pairs.  Otherwise,  finite  element  techniques  are  to 
be  introduced  for  the  deformable  elements  and  the  Jacobian  should 
be  constructed  accordingly.  Then  all  types  of  analyses  treated 
here  can  be  performed  in  a  routine  manner.  However, 
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for  average-shaped  mechanisms,  introduction  of  def ormability  through 
spring-dampers  is  recommended.  For  structural  problems,  when  the 

forcing  functions  depend  on  time  and  on  solution  variables,  the 

\ 

conventional  way  of  obtaining  global  solutions  through  modal  analysis 
and  Duhamel  Integrals  [54]  fails  and  the  techniques  employed  in  this 
development  may  be  employed.  For  structures  such  as  automobiles, 
where  finite  element  techniques  are  to  be  used  for  chassis  and  frames, 
in  conjunction  with  the  motion  of  rigid  body  mechanisms  in  the  sus¬ 
pension  systems,  sparse  matrix  techniques  and  stiff  integration 
(Gear)  algorithms  are  indispensible . 

Looking  back  into  the  details  of  the  DADS  computer  program,  it 
is  the  opinion  of  the  author  that  some  modifications  and  refine¬ 
ments,  particularly  in  the  field  of  numerical  integration,  are  nec¬ 
essary  to  make  it  more  efficient.  Among  other  refinements,  the 
following  are  under  active  consideration: 

(1)  Storage  of  the  values  of  only  those  solution  variables  and 
related  data  in  a  disk  in  transient  analysis,  that  are  to  be 
used  in  adjoint  analysis.  This  then  becomes  a  problem- 
dependent  affair. 

(2)  Allotment  of  only  one  disk  for  storing  the  solution  data  set 
of  the  adjoint  equations  corresponding  to  all  violated  con¬ 
straints  in  the  sensitivity  analysis. 


144 


APPENDIX 

MATHEMATICAL  RELATIONS  USED  IN  THE  SUBROUTINES 
RELATE,  DGDBZ ,  ADLDZ ,  AND  DLDFB 
FOR  THE  EXAMPLE  PROBLEMS 


In  the  following,  symbols  of  Chapter  II  and  V  and  data  given  in 
Chapter  VI  have  been  used. 

The  Slider-Crank  Mechanism 

Subroutine  RELATE: 

From  Fig.  2.5  and  the  initial  data  (Table  6.1),  one  obtains 


J3  =  3  m3 


b 


2 

3 


YF1(1)  =  b2 
YF2 (1)  -  b2  -  0.5 
Xl(3)  =  b3 


X2(2)  =  -b3 
X(3)  =  b3  cos  <j>3 


X(4)  =  2. DO  x  X(3) 
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Subroutine  DGDBZ : 

From  Eqs.  (6.6)  to  (6.9) 

DGDB(l)  =  DGDB(2)  =  DGDB (3)  =  DGDB (4)  =  0 

for  all  the  four  functional  constraints.  Only  for  the  cost  function 

DGDB (4)  =  1.0  and  DGDB (I)  =0,  I  =  1,2,3  . 

Also  for  the  fourth  constraint,  the  only  nonzero  component  for  DGDZM 
is  given  by 

DGDZM(9)  =1.0  . 

Subroutine  ADLDZ: 

From  Eq.  (6.6),  one  obtains  for  the  first  constraint  (since  as 
the  26-th  solution  variable) , 

DLDZ (26)  *  -Al* (1. DO  +  DSIGN (1.D0.A1) ) /10.D0 


where 


Al  -  - 


10. DO 


+  1.D0 


+  e  , 


e  being  the  e-active  constant.  All  other  compontnets  are  zero. 
D  after  decimal  point  indicates  double  precision.  The  function 
DSIGN(a^,a2) ,  a^,a2  being  the  arguments, is  defined  as  follows: 

DSIGN(a1,a2)  =  laj  sgn(a2)  , 


where 
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(  1 

for 

a2  >  0 

sgn(a2)  = 

for 

a2  =  0 

(-1 

for 

a2  <  0 

(6.7),  one 

obtains 

for 

the  second 

component  as 

DLDZ (26)  =  Al* (1.D0  +  DSIGN(1.D0,A1) ) /10.D0 

where 


A1  '10^0  -  1-D0  +  £ 


From  Eq.  (6.8),  one  obtains  for  the  third  constraint  the  only  nonzero 
component  as 


DLDZ (33)  =  2.D0*b1*Al*A2*(l.D0+  DSIGN(l.D0,A2))/b4 


where 


b  *  (Al)2 

A2  =  -  -  1.D0  +  e 

b/. 


From  Eq.  (6.9)  all  components  of  DLDZ  are  zero  for  the  fourth  con¬ 
straint. 


Subroutine  DLDFB: 

The  nonzero  components  of  DLDB,  DFDB,  DPHDB,  and  DZIDB  are  given 
by  the  following  relations.  For  the  third  constraint 


-f 


DLDB(l)  =  (A1)2*A2*(1.D0  +  DSIGN(1.D0 ,A2) ) /b^  , 

DLDB (4)  =  -A2*(1.D0  +  DSIGN(1 .DO, A2) ) *b'* (Al) 2/ (b^) 

where 


b  *  (Al)2 

A2  =  -  -  1.D0  +  e 

b4 

From  the  equations  of  motion,  one  obtains 

DFDB(15, 3)  =  ~(y3  +  y5)  sin  <f>3  +  (y^  +  y&)  cos  <f>3 

DFDB(21, 2)  =  F^  cos  <(>  +  Fb  sin  <J>4 

DFDB(3,2)  =  -F*  cos  4>  -  F*  sin  <j> 

A  XX  X 

From  the  equations  of  constraints,  one  obtains 


DPHDB(3, 3)  =  cos  <(>  =  DPHDB(5,3) 

DPHDB(4,3)  =  sin  <J>3  =  DPHDB(6,3) 


From  the  spring-damper  relations,  one  obtains 


DZIDB(3, 1) 


*  XI 


DZIDB(4, 1)  = 


*  YI 
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DZIDB(1, 2)  = 


[XI*(-sin  <j>^  +  sin  <p^)  +  YI*(cos  ^  -  cos  ^)] 


2  2  1/2 

[(xi)  +  (yd  r/z 


^1  ^ 

DZIDB (3 , 2)  =  b1  — ~ — —  *  (-sin  ^  +  sin  4>4) 


DZIDB(4,2)  = 


ai  -  ZP 

bx  - - -  *  (cos  ^  -  cos 


where 


XI  =  XG1  -  XG2,  YI  =  YG1  -  YG2 
XG1  =  XF1(1)  *  cos  <J>1  -  YF1(1)  *  sin  $  +  X(l) 

YG1  =  XF1(1)  *  sin  *  +  YF1(1)  *  cos  tp  +  Y(l) 

XG2  =  XF2(1)  *  cos  <j>4  -  YF2 (1)  *  sin  <j»  +  X(4) 

YG2  =  XF2 (1)  *  sin  $  +  YF2(1)  *  cos  <f>4  +  Y(4) 


c 


The  Plow-Share  Mechanism 
(with  re-entering  angle  of  0.0174533  radians) 

Subroutine  RELATE: 

From  the  initial  data  and  definition  of  the  design  parameters 
(Chapter  VI),  one  obtains  (see  Eqs.  (6.14)  and  (6.12), 


K1  =  23.25  x  IQ8  x 


b3  ^ 
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£0  =  bl  b3  +  °* 56604  *  10"2 
Subroutine  DGDBZ: 

In  this  case  the  only  nonzero  component  of  DGDB  is  given  by 
DGDB(4)  =  1.0 

for  the  cost  functional  only. 

All  the  components  of  DGDZM  are  zero.  ■ 

Subroutine  ADLDZ: 

The  nonzero  components  of  DLDZ  are  given  by  the  following.  From 
Eq.  (6.18),  one  obtains  for  the  first  constraint, 

DLDZ (12)  -  -A1  x  (1.D0  +  DSIGN(1.D0,A1) )  ,  for  t  >  0.6DO 

where 


A1  =  -<j>2  -  0.0174533D0  +  e 


From  Eq.  (6.19),  one  obtains  for  the  second  constraint 
DLDZ(57)  =  2. DO  *  K1  *  A1  *  A2  *  (1.D0  +  D  SIGN(l.D0,A2))/b4 


where 


A2  = 


1  2 
K  *  (Al) 


1.D0  +  e 
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Subroutine  DLDFB: 

The  nonzero  components  of  DLDB,  DFDB,  DPHDB,  and  DZIDB  are  given 
by  the  following  relations.  From  Eq.  (6.14) 


=  93.D08  * 


b3  ^ 


-69.75D08 


b3b 


4 

2 


-23.25D08 


b32b 


3 

2 


Then  for  the  second  constraint 


2  1 

DLDB(l)  =  A3  *  ((Al)  *  ~~  -  2. DO  *  K  *  Al  *  b,)/b.  ' 

3b^  3  4 

DLDB (2)  =  A3  *  (Al)2  * 


2  r\K  1 

DLDB (3)  =  A3  *  ((Al)  *  ~ -  -  2. DO  *  K  *  Al  *  bj/b, 

o  14 


DLDB (4)  =  -K1 


*  A3  * 


f 


where 
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A2 


1  2 
K  *  (Al) 


1. DO  +  e 


A3  =  A2  *  (l.DO  +  DSIGN(1.D0,A2) ) 


From  the  relations  corresponding  to  the  first-spring-damper  pair,  one 
obtains 


DZIDB(3,1)  =  DL1  * 


DZIDB(3,2)  =  DL1  * 


DZIDB(3,3)  =  DL1  * 


DZIDB(4,1)  =  DL2  * 


DZIDB(4, 2)  =  DL2  * 


DZIDB(4,3)  =  DL2  * 


3K 

3b 


ax; 

3b 


IE 

3b- 


3K 

3b, 


3K 

3b„ 


3K 

3b- 


DL3  *  b3 


DL3  *  b 


DL4  *  b3 


DL4  *  b 


where 
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DL1  = 


(£i  -  V 


*  XI 


DL2  = 


(£!  “  V 


*  YI 


DL3  = 


V1  *  XI 
K  *  — 

Z1 


DL4  = 


K1  *  T 

1 


XI  =  XG1  -  XG2  ,  YI  =  YG1  -  YG2 


XG1  =  XF1(1)  *  cos  <j>4  =  YF1(1)  *  sin  $  +  X 


YG1  =  XF1(1)  *  sin  $  +  YF1  (1)  *  cos  cj>4  +  Y^ 


XG2  =  XF2(1)  *  cos  -  YF2(1)  *  sin  <t>,  +  X, 

D  DO 


YG2=  XF2 (1)  *  sin  <f>6  +  YF2(1)  *  cos  <j>  +  Y& 
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